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ABSTRACT 


This research has dealt with the analysis of the in- 
depth response of phenolic-nylon and silicone elastomer 
ablative composites. The calculations of the total energy 
absorbed were performed by the coupling of the virgin 
plastic zone with the decomposition and the char zone. 

This was achieved by taking into account the energy ab- 
sorbed by the reacting gases that flow into the char zone 
as a result of the depolymerization of the virgin plastic 
ablator and/ in addition, by taking into account the 
energy absorbed associated with this depolymerization 
process. The bulk of the research was concentrated on 
phenolic-nylon ablative composite since it has shown the 
superior performance in comparison with other ablators in 
high heating rate environments. Frozen, equilibrium and 
non-equilibrium analyses were developed for this composite. 

A non -equilibrium or kinetic analysis was not developed 
for silicone elastomers, because of insufficient data on 
the complex solid-state reactions that are known to occur 
in the silicone-carbon system. However, for the silicone 
elastomers the in-depth response was modeled by using 
equilibrium analysis which demonstrated the flexibility 
and generality of the computer analyses developed. 

V 

In this research we were able to identify the important 
energy absorbing mechanisms in an ablator. In the plastic 
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zone the predominant energy absorbing mechanism is due to 
the decomposition of the plastic itself. In the char the 
predominant mechanisms are transpiration cooling , endother- 
mic chemical reactions, and sublimation of carbon. In 
conjunction with the equilibrium and frozen flow analyses 
it was established that chemical reactions become important 
energy absorbing mechanisms above 2000°F, while sublimation 
does not become an important contributing factor unless the 
temperature is above about 4700°F. 

For the phenolic-nylon ablator it was determined that 
a non-equilibrium flow analysis, employing detailed kinetic 
equations of the pyrolysis gases and a detailed analysis of 
the depolymerization of the virgin plastic composite, was 
necessary to accurately describe the in-depth response of 
the ablative composite. 

In the non-equilibrium flow analysis the important 

chemical reactions and kinetic data for a temperature range 
o o 

of 500 F to approximately 6000 F were determined by 
screening over one hundred chemical reactions. For each of 
the fifteen reactions finally selected to model the chemical 
behavior of the pyrolysis gases an equilibrium constant was 
computed to provide a means of calculating the reverse 
reaction rate constant and therefore, make the reactions 
reversible and thermodynamically consistent. 

The decomposition/char zone temperature interface for 
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all the cases analyzed was always higher than 2000°F in 
the non-equilibrium flow analysis, and the pyrolysis 
gases that transpired into the char zone reacted very 
quickly. As a consequence of the very rapid reaction rate 
in the non-equilibrium analysis it was necessary to use 
very small step-sizes in the Runge-Kutta numerical solution 
in order to maintain the stability of the solution. This 
resulted in excessive amounts of computer time which 
limited the number of cases analyzed. The numerical diffi- 
culties experienced with the non-equilibrium analysis 
were due to a phenomenon called stiffness. Because of this 
the minimum total mass flux analyzed for the non-equilibrium 
flow analysis was 0.35 lb/ft 2 -sec. However, no such res- 
trictions were necessary for equilibrium or frozen flow. 

Parameter studies with pressure and char surface 
recession velocities were conducted to determine the effect 
that these variables had on the total energy absorbed by 
the ablator. It was found that pressure had very little 
effect on the energy absorption below 4000°F, both for 
equilibrium and non-equilibrium analyses. Surface recession 
velocity had an effect on the decomposition/char zone 
temperature interface, with the interface temperature 
increasing with increasing recession velocity. 

X 

Although frozen flow analysis proved to be a very poor 
approximation to the total amount of energy absorbed. 
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equilibrium analysis, using free energy minimization, 
proved to be a reasonable approximation to non-equilibrium. 
Although it suffers from its inherent simplifying assumption 
that the gases are always in chemical equilibrium. 

It was found that the silicone elastomer ablator 
showed the same trends observed for the phenolic-nylon 
ablator for the several parameter studies considered. 
However, it was found that the phenolic-nylon ablator is 
the more efficient of the two because it absorbs more 
energy per pound of material. 
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CHAPTER I 


INTRODUCTION 

This ; research describes the decomposition in-depth of 
ablative composites along with the transport phenomena of 
pyrolysis gases which result from the decomposition of these 
plastics as they flow through the porous char of char- 
forming ablators. In particular, the pyrolysis products are 
those formed by the thermal degradation of nylon-phenolic 
resin and silicone elastomer composites. The nature and 
extent of chemical reactions among the pyrolysis products 
and the char, along with the energy absorbed by the combined 
pyrolysis and char zone, are given major emphasis in this 
research. Likewise, the determination of the important che- 
mical reactions with thermodynamically consistent kinetic 
data are necessary in developing a realistic analysis for 
predicting the thermal performance of ablative heat shields , 
and they have been obtained in this research. 

The Design of Thermal Protection Systems 

It is a well known fact that vehicles reentering 
planetary atmospheres need to dissipate enormous amounts 
of energy. This energy dissipation is achieved by reducing 
the speed of the vehicle by either the use of rockets or by 
taking advantage of the frictional drag of the atmosphere 
to decelerate the vehicle. The former solution is not an 

i 
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acceptable design practice since it increases the fuel 
requirements for a mission. This presents, however, some 
severe thermal and structural requirements for reentry 
vdiicles (1, 2) since all of the energy the vehicle posse- 
sses at reentry, kinetic as well as potential, must be con- 
verted into heat. To illustrate that this energy can indeed 
be great. Figure 1-1 compares the energy per unit mass as a 

function of entry speed V , in kilometers per second, for 

E 

several heat shield materials (3) . To dissipate this great 
amount of energy several methods of solutions are possible 
(4, 7, 8). These are: (1) Heat Sink, (2) Radiation Cooling, 

(3) Transpiration Cooling and (4) Ablative Cooling. 

Heat Sink : This method consists of providing enough 

mass with high enough heat capacity to safely absorb the 
heat input to the vehicle (12) . 

Radiation Cooling : This method of cooling is achieved 

by using highly reflective surfaces such that most of the 
heat incident to the body, by either convection and/or 
radiation, is reflected back rather than absorbed by the 
vehicle (13) . 

Transpiration Cooling : This technique injects fluid 

into the boundary layer through openings in the body. This 
injection of mass produces a substantial heat blockage that 
reduces the net heat transfer to the body by virtue of the 
thickening of the boundary layer. 

Ablative Cooling : This form of thermal protection is 
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FIGURE l-l COMPARISON OF THE ENERGY PER UNIT MASS 
AS OF FUNCTION OF ENTRY SPEED , (5). 
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achieved by the sacrificial loss of the material covering 
the body. The energy absorption is achieved by a change of 
phase. This phase change can be melting, vaporizing, and/ 
or subliming. 

Of the four techniques, the first three suffer from 
excessive weight penalties to the vehicle. The feasibility 
of ablation heat shields for satellite reentry was esta- 
blished in the late 50's and early 60's (2). Although 
absorption of heat by phase change is the distinguishing 
feature of ablation, energy dissipation by radiation, con- 
duction, convection, transpiration, cooling and chemical 
reactions is likewise achieved (8, 9). Georgiev, et. al. 

(5) has a discussion of the various factors affecting the 
choice of a heat sink, ablation, or radiation shield to 
use for thermal protection systems. 

In addition, Lees (11) , presents an example which illus- 
trates, that for a good thermal conductor such as copper, 
the peak heat transfer rate should be limited to about 1,000 
BTU/ft -sec, for 30 seconds for a maximum allowable tempera- 
ture of 1500°F. This clearly shows the limitations which 
are imposed on vehicles made of non-melting solids. Scala 
(14) compares the effectiveness of transpiration cooling with 
ablating cooling and concludes that ablation is more effec- 
tive than transpiration cooling on a coolant mass consump- 
tion basis. Later, Beusman and Weisman (15) showed that 
transpiration cooling would be more effective if water were 
to be used as the coolant. 
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Aerodynamic Heating of Blunt Bodies vs . 

Slender Bodies 

It was said that for a body entering the atmosphere, 
which is braked by aerodynamic drag, all of the energy 
possessed by the vehicle is converted into heat. Fortunately 
however, only a fraction of the total kinetic energy need be 
transferred to the vehicle as, heat. The approximate amount 
can be determined from the following analysis presented by 
Allen (3). He simplified the motion equation to consider on- 
ly the most important term, the drag term. Then, 


mdV 

“Ht 




( 1 - 1 ) 


Where m is the vehicle's mass, V E the entry speed, t the 
time, p the air density, A the reference area and C D the 
drag coefficient based on A. 

The rate of heat transfer to the vehicle was shown to 

be: 
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c H p V A 


( 1 - 2 ) 


where q is the heat input and C R is the heat transfer 
coefficient. 

By combining these two equations it was shown that the 
heat input for the entire braking process was given by the 
following expression: 



That is, the amount of heat given to the body is determined 



6 


by the ratio, which he called, the "energy ratio", of the 
heat transfer coefficient, C H , to the drag coefficient, C D . 
Using Reynolds analogy and assuming that only convective 
heating occurs, it was shown that: 


n 



c 


F 


2 (C Dp + C p ) 


(1-4) 


where C F is the skin friction coefficient and C Dp is the 
pressure drag coefficient. It should be noted that blunt 
bodies have much greater pressure drag coefficient than 
slender bodies. Hence the fraction of heat transferred to 
the body is less than for slender bodies. This interesting 
relation formed the basis for using blunt bodies as reentry 
vehicles under high heating load conditions and for long 
time reentry (6) . "Depending on the relative amount of 
laminar and turbulent flow a slender vehicle designed with 
a length to diameter ratio of 3.0, for example, would have 
between about 6 to 16 times as severe a heating problem as 
a blunt vehicle designed to produce no lift", (6). 


Thermal Protection System for Manned Reentry 

It was previously stated that heating for a reentry 
vehicle is produced by conversion of the kinetic energy of 
the moving body into thermal energy when it decelerates. 
For an efficient thermal protection system most of the in- 
cident heat must be disposed at the outer surface. Any of 
the four possible methods previously mentioned can be used 
effectively. However, it has been demonstrated that abla- 
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tion cooling is best suitable for manned reentry missions 
(3,8). 


Ablative Materials : Ablative materials are generally 

made of a ceramic, a plastic, or reinforced system formed 
by a combination of plastic and an inorganic fiber. In ge- 
neral plastics and related composites have been widely used 
(16) . Ideally, these ablative materials must possess low 
thermal conductivity, high heat capacity and large heats of 
degradation (8) . This is needed to effectively protect the 
vehicle from the high heating environment encountered during 
reentry. 

Success has been achieved by employing composites of 
nylon, phenolic resin, silicon elastomers and others. A 
partial list of some of these composites is presented in 
Table 1-1. These composites fall into broad categories, 
non-charring (thermoplastic) and charring (thermosetting). A 
non-charring ablator is one which vaporizes and decomposes 
into gases with little or no residue remaining on the abla- 
tive surface (17) . Teflon is one such non-charring com- 
pound. The charring ablator vaporizes, but it decomposes 
into low molecular weight pyrolysis gases and a carbonaceous 
residue. The gases are heated as they travel through the 
char layer with chemical reactions occurring. These gases 
when injected into the boundary layer effectively block part 
of the heat input to the vehicle surface by thickening the 
boundary layer. Phenolic-nylon is an example of a charring 



Table 1-1. List of Materials Tested for Heat Protection of Reentry Vehicles ( 14 ) 
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ablator. "In general charring ablators, subliming ablators 
and melting ablators provide the most effective thermal 
tection system" (18) . Of these, the charring ablator nor- 
mally provides the most efficient thermal protection shield 
for a wide variety of applications including manned reentry 
vehicles (18,19). For this reason charring ablators have 
been studied. 

/ 

Ablation ! " Ablation is an orderly heat and mass 
transfer process in which large amounts of thermal energy 
are expended by sacrificial loss of surface region material" 
(19) . To date certain physical and chemical aspects of the 
process are well known and understood. A description of 
those physical and chemical processes for a charring abla- 
tor follows. 

Charring Ablator; Physical Process : on high speed 

planetary reentry, the high temperature generated by vis- 
cous forces causes the plastic ablator to undergo a number 
of physical and chemical changes with the formation of a 
carbonaceous residue. A cross section of the ablative com- 
posite with the accompanied flow field is illustrated in 
Figure 1-2. 

The heat input from the surrounding flow field is ab- 
sorbed, dissipated, blocked, and conducted into the plas- 
tic material substrate (19). Initially, the heat flow 
in penetrates at a low rate because of the low thermal 
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convective heat transfer to the char surface 


Figure 1-2. Schematic Diagram of the Various Zones 
Developed During Reentry of a Capsule 
Protected by a Char Forming Ablative 
Heat Shield. 
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conductivity of the ablator. This causes a rapid 
temperature rise on the surface with subsequent phase 
change and thermal degradation of the material into 
pyrolysis gases and a porous carbon char. During this 
phase change of the virgin plastic composite into 
bases and char an enormous amount of heat is absorbed. 
This prevents the interior of the vehicle from getting 
excessively hot. The pyrolyzed gases flow through the 
char to the surface and in the process absorb heat due 
to their own heat capacity (transpiration cooling) and in 
addition, undergo, mostly endothermic chemical reactions. 
The pyrolyzed gases are injected into the boundary 
layer producing a blowing effect and thus, reducing the 
net heat transfer to the surface of the vehicle. The 
char produced from the decomposition of the plastic 
material and the further cracking of the pyrolysis 
gases, has a double protective effect. It acts as an 
insulator for the virgin material by effectively reducing 
the rate of heat conducted to it. It also acts as a 
radiation shield by radiating back part of the heat 
incident to the surface. Thus the char serves to control 
the surface temperature and greatly restrict the flow of 
heat into the substrate interior (19) . 

Summary 

This chapter has illustrated the severe heating 
problems encountered by vehicles entering planetary 



atmospheres. It was shown that for high heating rates 
and long reentry time blunt body vehicles would absorb 
less heat than slender body vehicles. Of the several 
methods proposed to overcome the thermal barrier, 
ablation cooling was shown to be superioTf to the others 
for manned vehicle applications. In addition, it was 
shown that char forming ablators would normally provide 
better thermal protection for a wide variety of applica- 
tions including manned reentry vehicles. Finally, it 
was illustrated that most of the heat transferred into 
the vehicle was absorbed at the decomposition zone and 
by the gases flowing through the char zone. For this 
reason, it is important that any theoretical analysis 
designed to predict the thermal response of ablative 
composites have an accurate description of the de- 
composition zone and char zone since they are the im- 
portant heat absorbing regions. In the subsequent 
chapter several analyses proposed to describe the 
decomposition zone and char zone are reviewed in the 
light of present knowledge of the ablation process. 
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CHAPTER IX 

A REVIEW OF ANALYSES DESCRIBING DECOMPOSITION IN- 
DEPTH FOR THE ABLATIVE PROCESS OF CHAR 
FORMING ABLATORS 

Research and development of ablative heat shield compo- 
sites for space vehicles can be grouped into two broad 
categories. The first involves the detailed investigation 
of the physical and chemical processes which occur during 
ablation. In these studies particular attention has been 
given to the experimental investigation of plastic decompo- 
sition chemistry (1-4) , valid analytical descriptions of the 
pyrolysis zone (5-10), and the flow of pyrolysis gases 
through porous media (11-22) . In addition, studies have 
been made of the effects that char oxidation (23-32) , 
thermal property variations (33-39) , porosity and permeabi- 
lity of the char (34-37, 17), composition of the pyrolysis 
gases entering the char zone (40-43 ) » environmental condi- 
tions (43-47), char spallation (48-50), boundary layer 
interaction (51-58) , radiation (59-69) and others (70-92) , 
have on the accurate prediction of the thermal performance 
of ablative composites. The second category covers the 
analysis of the transient response of the combined heat and 
mass transfer mechanism which occurs between the heat shield 
and the flow field (7, 51, 53, 58, 84, 111). Research in 
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both areas is essential in developing more effective ther- 
mal protective system for reentry vehicles. The former 
category improves the accuracy of the transient response 
computation by supplying better experimental and theoreti- 
cal understanding of individual physico-chemical processes 
which occur during ablation. The latter allows for a more 
economic and effective thermal design of the ablative heat 
shield. 

As mentioned in Chapter I, this research describes 
the decomposition in-depth of ablative composites along 
with the transport phenomena of pyrolysis gases which re- 
sult from the decomposition of these plastics as they flow 
through the porous char of char-forming ablators. In 
particular, the pyrolysis products are those formed by the 
thermal degradation of nylon-phenolic resin and silicone 
elastomer composites. The nature and extent of chemical 
reactions among the pyrolysis products and the char, along 
with the energy absorbed by the combined pyrolysis and char 
zone, are given major emphasis in this research. Likewise, 
the determination of the important chemical reactions with 
thermodynamically consistent kinetic data are necessary in 
developing a realistic analysis for predicting the thermal 
performance of ablative heat shields, and they have been 
obtained in this research. 

In one phase of this research April (19) made a state 
of the art study of the flow of pyrolysis gases in the char 
zone and reviewed the important works related to this area. 
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In this chapter the emphasis is placed on the quasi-steady 
analysis for describing decomposition in-depth of virgin 
material composites. In addition, some of the typical 
transient response analyses are briefly discussed. In 
particular, the transient analyses of Swann, et. al . (58) , 
Clark (111) , Kendall, et . al. (51,53) , and Kratch, et . al . 
(83), will be reviewed. Finally, a summary of previous 
research on flow in the char zone are presented. The dis- 
cussion of analyses used for the description of the decompo- 
sition process in virgin plastic ablative composites follows 
below. 

Decomposition In-Depth of Virgin Plastic Composites 

Several methods are available for the treatment of the 
thermal decomposition process of plastic composites. They 
differ primarily in whether the chemical decomposition of 
the plastic is assumed to occur in a single reaction plane 
(5, 6) at a fixed temperature, or whether a spacially 
continuous decomposition in-depth is used (7, 8, 9). The 
analysis in this research assumes the decomposition of the 
plastic ablator to occur spacially, and the temperature 
limits within the reaction zone is determined by the kine- 
tics of the decomposition of the virgin material. The 
frequency factor, activation energy and the order of 
equations can be obtained from thermogravimetric analysis 
for particular materials of interest (1-4) . In Figure 2-1 
a typical thermogravimetric curve is shown for a low density 
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Figure 2-1. Typical Thermogravimetric Curve Relating the Mass 
Loss as a Function of Temperature for a Phenolic- 
Nylon Resin (2). 
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nylon-phenolic resin composite (2) . 

It is possible to calculate from this data the rate of 
pyrolysis gas generation W^. The rate of mass loss of a 
material is affected by the rate at which the temperature is 
increased (deg/sec) ; that is, the heating rate. The 
inaccuracies resulting from the effect that the heating rate 
has on the pyrolysis gas generation rate, W , can be elimina 
ted by the use of a kinetic expression of tne form of 
Equation (2-1) below (1). This equation describes the rate 
of change of density, with time and is independent of the 
heating rate. 


-1 dpj 
p i,o dt 


<»± - Ps)/ p i,o> ni *e‘ (E i /RT) (2-1) 


One approach to modeling the decomposition of plastic 
composites has been presented by Kondo and co-workers (3). 
Their analysis assumes decomposition to occur spacially 
rather than on a single plane. The authors state that this 
method of analysis should be used when the plastic material 
is made of different composites. They compare their analy- 
tical solutions with experimental data and report good 
correlation. 

Another study of the decomposition process for plastic 
composites has been presented by Stroud (5, 6). In this 
investigation, the effects that kinetics, heat of pyrolysis, 
and surface recession velocity has on the reaction zone 
thickness, and the temperature range at which decomposition 
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of plastic occurs is examined. Kondo and co-workers (8) / as 
mentioned, noticed also that the reaction zone thickness is 
affected by different parameters. In an attempt to simplify 
a very complex process Stroud (5, 6) studies, in addition, 
the assumption that the decomposition of plastic composites 
could occur at a single plane. His analysis shows that this 
simplification is usually valid for one reacting species and 
is considered to be a good engineering approximation. 

Another analysis which considers decomposition reactions 
to occur over a finite region is that of Scala and Gilbert 
(9) . Their analysis differs from that of Kondo et . al . (8) 
and Stroud (5, 6) in that the reaction zone thickness is 
assumed to be a constant. This is a drawback since this 
does not allow for the investigation of the effects of kine- 
tics and heating rate on reaction zone thickness. Scala 
and Gilbert (9) discuss, in addition, some of the difficul- 
ties encountered by earlier investigators in describing 
the decomposition process accurately. They attributed the 
chief difficulty to an experimentally observed phenomena 
which showed that the mass generation rate, W^, varied 
widely at high heating rates. Better experimental techni- 
ques and better correlation of experimental data (1, 2) 
seemed to have diminished some of the difficulties encoun- 
tered by earlier investigators. An equation of the form of 
(2-1) was suggested by Friedman (9) to eliminate the 
uncertainties of the mass generation term caused by high 
heating rates. 
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Since decomposition in-depth is an important ingredient 
in the analysis of the ablation process, the work of Stroud 
(5, 6) and Kondo and co-workers (8) are reviewed in this 
chapter in detail. The work of Stroud (5, 6) is important 
because it determines and examines the important parameters 
which affect decomposition and reaction zone thickness. It 
is significant, in addition, because it attempts to simpli- 
fy a very complex process, that of the decomposition of 
plastic composites, with the assumption that decomposition 
occurs at a single plane, find shows under which conditions 
this simplification holds. 

The work of Kondo and co-workers (8) is also important 
in that they establish the necessity of considering decom- 
position to occur over a finite region and they examine, as 
in the case of Stroud (5, 6), the important variables that 
affect reaction zone thickness. The work of Stroud (5, 6) 
is presented in the following section. 

Study of the Chemical Reaction Zone in Charring Ablators 
and the Reaction Plane Approximation During Thermal De-~ 
gradation by Stroud . 

The work (5, 6) was a systematic approach of studying 
and establishing the important variables that affect the 
theoretical treatment of the thermal decomposition of 
plastic ablative composites. In addition, the validity of 
the approximation of describing decomposition of the plastic 
as taking place in a single reaction plane was studied. 

This was found to be a good engineering approximation for 
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one reacting species. The analysis with its assumptions 
is presented, and the results are summarized below. 

The energy flow in the ablative heat shield on a large 
body entering the atmosphere is approximately one-dimen- 
sional, and therefore the energy equation for a chemically 
reacting solid was approximated by Stroud (5, 6) as: 


c 

Z 3 Pi 


LL± Hi + p, C D , 3T 

i-l 3tT X * B . 1 p i/S 


3T - 


3W, 

-ri H i.g 

3y 


W.c_ 3T = 3 (k_ 3 T) 

p i,g Sy W 8 3y 


( 2 - 2 ) 


The above equation carries the assumption that the kinetic 
energy term can be neglected since it is three orders of 
magnitude smaller than the combined effects of the chemical 
and thermal energy terms. In addition the pressure was 
assumed to remain constant throughout the material because 
of the high permeability of chars. The author reasoned 
that "in practice, materials which form chars of low per- 
meability and low strength would not be used due to the 
danger of internal pressure blowing off periodically and 
thus reducing the efficiency of the material" (5). 

From the continuity equation 



V • pu 


(2-3) 
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where u is the velocity vector for the pyrolysis gases, 
and, therefore, for one-dimensional flow it becomes: 

3p = - 3W (2-4) 

<Tt Sy 

Substitution of Equation (2-4) into (2-2) results in: 


c 


l 

i*=l 



3T = 

<nr 


3 (k 

Ty 


a 


3 T) 

3y 






(2-5) 


Notice that the number of gas components is numerically 
equal to the number of solid components in the ablator. 

This stems from the assumption made (5, 6) that each solid 
component degrades to a single gas species. 

By assuming a quasi-steady state, that is, one in 
which the total mass flux of the system is constant, Stroud 
was able to eliminate the time dependence of the equation 
by using the following transformation: 


y = y - vt 


where v is the surface recession velocity. Therefore, after 
the transformation is performed Equation (2-5) becomes: 

k 3 + ( f< c P ( «i as - v a Pi 4 Hjl + 
dy2 i=l i,g dy 


dT ) 


v p i c p, s 

1,8 dy 


( 2 - 6 ) 
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It should be noted that the assumption of constant 
thermal conductivity is valid only over a short tempera- 
ture range, since the term dk/dT can become appreciable 
at higher temperatures. In our research the restric- 
tion of constant thermal conductivity was not made. 

Stroud (5, 6) further simplified the energy equa- 
tion by assuming that the specific heats of all 
constituents were the same, and from the quasi-steady 
assumption: 


W = vp Q - vE 

i=l 


(2-7) 


where P Q is the density of the virgin composite before 
degradation. Equation (2-6) was put into the form: 
c 


k a d*T Z 


I (-vdp j AH, + v C n dT)=* o 

, l-z -pr- ° P i,s 5 ^ 

dy2 dy y 


(2-8) 


By integrating Equation (2-8) once, evaluating the constant 
of integration and letting Ti=y Cp Po/^s ' Equation (2-8) 

3 . , 8 

was transformed to: 


c 

dT + E 
dr) i=l 


AH i p i,o (1- p i ) 


L p i. 


p i,o 


+ T-T q * 0 


(2-9) 


To obtain Pi as a function of the temperature T, and the 
dimensionless distance, n# an Arrhenius type relation of 
the form shown below was used: 


dpi 

cTF" 


- Pi n A^" E i/RT 


(2-10) 
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where is a dimensionless density of the form: 

! i 

Pi « P 1 (2-11) 

p i,o . 

which expresses the ratio of the density of the composite 
i at an instant of time with respect to the original density 
Q , of the virgin material. A^ is the frequency factor 
of the reaction and n, the order of the reaction. 

It is worthwhile to note at this point, that Stroud 
(5, 6) assumed that the degradation of each composite i can 
be described by a single kinetic expression. However, it 
should be noted that this is not true in most cases. For 
phenolic-nylon and for other composites as well, more them 

i 

one Arrhenius expression is necessary to describe the de- 
composition kinetics (1, 2) . 

From the continuity equation we have: 

dp i = - V d Pi (2-12) 

Sy"" 

and using the previous definition of n we have that Equation 
(2-12) becomes: 


d p i 

“dt 



v 2 dp ± 

XT 


(2-13) 


which when substituted into Equation (2-10) gives: 
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- n 

P, A.e 


-Ej/RT 


I# s 


P Q v' 


d Pi 


(2-14) 
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which expresses the change of the dimensionless density 
Pj_ / with the dimensionless distance, n* 

Two expressions were obtained by integrating Equation 
(2-14). One for n ^ 1 and was shown to be: 


p i = 


(n-1) 


^s 

C *i,s P ° v2 


- (Ei/RT) 


/RT) 1 l-n ± 
•dn + 1 J 


(2-15) 


h=o / 

The other was for a reaction of order one; i.e., n=l, and 
was shown to bes 


Pi = exp 


's 


' p i, 




s 



(2-16) 


At this point in his development Stroud (6) simplified his 
analysis by considering only that one degradation reaction 
was taking place. For this case, the following transfor- 
mation was found to be convenient.' 


0 *= RT (2-17) 

IT 

Applying this transformation to Equation (2-9) a 
dimensionless form of the energy equation was obtained. 


d0 + RAHp 1 (l-p) + 0 RT o = 0 (2-18) 

where p-^Q, is the initial density of the only reactable 
species present. 

Equations (2-15) , or (2-16) and Equation (2-18) were 
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solved simultaneously by numerical techniques to determine 
the value of temperature as function of distance in the 
ablation material. The numerical scheme used, as reported 
by Stroud (6), consisted of first holding the value of the 
temperature constant in the density equation and integrating 
the equation over a small increment of distance. Then, 
the density was held constant at its calculated value in 
Equation (2-lj3) while this differential equation was solved 
over a small spaoial increment. 

I | 

Stroud (6) solved the energy equation for a number of 
different conditions by varying such parameters as the 
frequency factor, the activation energy, the heat of pyro- 
lysis, the surface recession rate and the order of the 
reaction, and studied the effect that these parameters had 
on such things as the temperature profile reaction zone 
thickness and the reaction plane approximation. In Figure 
(2-2) the local density through the reaction zone is shown 
as a function of distance for two different values of the 
frequency factor. The curves shown were obtained for first 
order reactions. These curves were arbitrarily translated 
so that they would coincide when 50% of the degradable 
plastic had decomposed. 

The reaction plane approximation is also shown. The 
char thickness was determined by considering that the char 
recession rate and the linear velocity of the reaction 
zone were the same. The thickness of the reaction zone was 
arbitrarily defined as the distance between the point at 
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Figure 2-2. Variation of the Density Profile 

of the Virgin Plastic Through the 
Polymer Reaction Zone with Frequency 
Factor, as Reported by Stroud 0>). 
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which the local density reached a value of 98 percent and 
the point at which it reached 2 percent of the original 
density of the reactable species. In the example shown by 
Stroud (6) 80 percent of the material was reactable. The 
two values of the frequency factor used in Figure 2-2 were 
said to be typical of a number of ablation materials. An 
increase in the frequency factor from 10® ^ec“* to 10 *2 sec”* 
produces a decrease in the reaction zone thickness of about 
50 percent, from 2.1 mm to 1.6 mm. This seems reasonable 
to expect since a higher frequency factor implies a greater 
degradation rate and, therefore, a shorter time for degra- 
dation with a corresponding smaller reaction zone thickness. 
The two values of the frequency factor shown are said to be 
typical of a number of ablative materials. 

Figure 2-3 exhibits the effects of surface recession 
velocity on reaction zone thickness for two values of the 
heat of pyrolysis. The curves were calculated for first 
order reactions and a frequency factor of 10*2 sec _ *. a 
recession velocity of 25 micrometers per second (ym/s) was 
said to be representative of those obtained for charring 
ablators presently known. The curves of Figure 2-3 show 
that surface recession velocity has a greater impact than 
heat of pyrolysis on reaction zone thickness. The reason 
why the reaction zone thickness decreases with increasing 
surface recession velocity is that the greater the surface 
recession velocity, the greater the temperature gradient 
across the reaction zone. This steep temperature gradient 
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causes the plastic to degrade rapidly thus reducing the time 
of reaction and hence, the thickness of the reaction zone. 

It appears that if the reaction zone thickness decreases 
with increasing surface recession velocity, the reaction 
plane approximation is asymptotically approached as the 
surface recession velocity increases. And therefore, 
the reaction plane approximation becomes more accurate the 
higher the surface recession velocity. An increase in the 
heat of pyrolysis produces a decrease in the thickness of 
the reaction zone. However, the size of the change of the 
reaction zone thickness is not as great as with the surface 
recession velocity. 

Figure 2-4 illustrates the effect of activation energy. 
It demonstrates that the reaction zone thickness is practi- 
cally independent of activation energy. "These results 
lead to an important conclusion regarding the effects of 
the frequency factor on ablative performance. The thickness 
of the reaction zone is controlled almost entirely by the 
frequency factor, whereas the pyrolysis temperature is 
strongly influenced by activation energy" (6) . Activation 
energy and frequency factors have been shown to affect 
pyrolysis in different ways, hence a unique value of each 
must be found if accurate results are required. 

Stroud (6) proceeded to define, what he called, the 
"median temperature of reaction” as that temperature at 
which 50 percent of the reactable material has been degra- 
ded. By using the median temperature definition, to get 
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that point in space where the degradation reactions were 
supposed to take place in a plane, he was able to obtain 
correlations for the mass generation rate for reactions of 
order one-half, one and two. These are shown below: 
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for second order reaction: 
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where the dimensionless median temperature 1/0 is obtained 
by using the equations shown below. . These equations were 
developed by Stroud (6) and are: for a half order reaction: 




= -0.7 + 0.955 In A - 2.2 AH 


for a first-order reaction: 

1 = - 1.50 + 0.955 In A - 3.6 AH 


m 


and for a second order reaction: 


(2-22A) 


(2-22B) 


= - 1.70 + 0.945 In A - 1.8 AH 


(2-22C) 
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Equations (2-19) , (2-20) and (2-21) hold when 

< 

O' < 4H < 0.5 
iO 6 < A < 10 15 

Equation (2-19) , (2-20) and (2-21) give the mass flow 
rates of pyrolysis gases for the reaction plane approxima- 
tion for one reactable species. These flow rates locate 
the reaction plane at the median point in ^the reaction zone 
as was shown in Figure 2-2. Thus, it was shown that the 
rate of pyrolysis is a function of thermophysical data 
which can be readily obtained from thermogravimetric ana- 
lysis. The correlations derived by Stroud (6) are only 
valid for one reactable species, and of course for high 
surface recession velocities. The author did not present 
a study, nor did he speculate on what possible effects 
more than one reactable species would have in his analysis 
of the reaction plane approximation. However, he did show 
a comparison between the reaction plane approximation and 
the more complicated in-depth analysis for the case of one 
reactable species. This is illustrated in Figure 2-5. For 
each frequency factor shown, the temperature distribution 
obtained with the reaction plane approximation is close to 
that obtained with the more complicated in-depth analysis. 
This agreement, as stated by the author, illustrates the 
accuracy that is possible to achieve using the reaction 
plane approximation when quasi-steady conditions exist. 

Several conclusions were drawn by Stroud (6) from 
studying and solving the equations for a wide range of 




Figure 2-5. Comparison of Temperature Profiles 
Calculated by Reaction In Depth 
Analysis and Reaction Plane 
Approximation (5, 6). 
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quasi-steady state conditions. These were: 1) The tempera- 
ture distribution in the ablation material is strongly 
dependent on the heat of pyrolysis; 2) Since activation 
energy and frequency factors have different effects of 
pyrolysis, it is necessary to obtain unique and accurate 
values for each of these parameters; 3) The total mass 
flow of pyrolysis and the median temperature can be 
correlated on the basis of the Arrhenius formula; 4) The 
reaction plane approximation, which incorporates the 
Arrhenius pyrolysis temperature , gave satisfactory agree- 
ment with reaction in-depth analyses and with experiments 
for one reactable species. 

The contributions of Stroud's studies (5, 6) have been 
to define those conditions under which the reaction plane 
approximation can be used for one degrading species 
without loss of accuracy. In addition, his studies have 
shown the critical effect that decomposition kinetics have 
on the accurate description of the ablation process. This 
indicates the necessity of obtaining unique and accurate 
values for the frequency factor and activation energy of 
the degradation process. In the following section the work 
of Kondo and co-workers (8) is discussed. 

An Analysis of Steady State Ablation of Charring Materials 
by Kondo, E*ujiwara and Matsomuto 

Physical Model : The work of Kondo, et. al. (8) is a 

one-dimensional, steady-state, ablation analysis where the 
thermal decomposition reaction in the charring material is 
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treated separately from the aerothermochemistry of the 
flow adjacent to the surface and the gas injection from 
the surface of the solid material. The study "is restric- 
ted to the analysis of the pyrolysis reaction in the 
charring material". (8). The material is considered to 
consist of three chemically different zones, namely, 
virgin material (solid) , pyrolysis zone (solid and gas) and 
the char layer (solid and gas) . In Figure^ 2-6 these three 
zones are illustrated. The pyrolysis zone regression 
velocity is assumed to be equal to the surface recession 
of the char. The recession of the char layer is assumed 
to be caused by three mechanisms: (1) Char layer being 

compressed by stagnation high pressure, (2) Char being 
consumed by the reaction of the pyrolysis gases with the 
carbonaceous char, and (3) by heterogenous reactions with 
the gas on the surface. 

Basic Assumptions of Hondo's et. al. Analysis : Before 
introducing the basic equations of the analysis we shall 
summarize some of their basic important assumptions: 

1) One -dimensional analysis. That is, the heat 
flow and gas flow were assumed normal to the ablating 
surface. 

2) Quasi-steady analysis. Thus the recession veloci- 
ty, v, is constant. 

3) The pyrolysis reaction is a simple chemical 
reaction which changed the plastic material into gases plus 
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Figure 2-6. A Schematic’ Diagram of the Three 
Chemically Distinct Zones Used by 
Kondo et, al. (8) in Their Analysis 
of SteacTy STate Ablation of 
Charring Materials. 
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char. 

4) Thermal equilibrium between the gas arid solid was 
assumed. 

5) The assumption of constant physical properties 
across the reaction zone was used. 

6) The viscocity of the gas was neglected in the 

energy equation . j 

7) The thickness of the ablating material was consi- 
dered to be large in comparison with the width of the char 
layer or pyrolysis layer. Therefore, the virgin material 
was assumed to be located at the infinite distance. 

It is pertinent to comment at this point on some of 
the above assumptions since they relate directly to the 
analysis proposed in our research. The assumption of 
thermal equilibrium is acceptable. Clark (16) found in his 
experimental studies with simulated chars that this was a 
fair approximation of what actually occurs during the 
ablation process. April (19) in his experimental work 
found in addition, that there was less than 200° F between 
simulated pyrolysis gases entering a char and the carbo- 
naceous matrix. This temperature difference was measured' 
approximately 1/4" away from the back surface of the char 
before the gases entered it. The assumption of thermal 
equilibrium allows in addition a simplification of the model 
by not having to consider the energy transfer between the 
gas and the solid matrix. The decomposition of the virgin 
material is assumed to be described by a simple Arrhenius 
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formula. As was mentioned before, this approach tends to 
make the kinetics independent of the heating rate at which 
they were obtained. This approach has been used also by 
Stroud (5, 6), Scala and Gilbert (9) and others (10, 84, 

85) . in this research the Arrhenius expression is also 
used. However, a more complicated model of seven kinetic 
reactions is used. Kondo and co-workers (8^) neglect the 
viscosity of the gases and, therefore, assumed that viscuous 
effects on the energy equation are negligible. In Chapter 
HI this assumption is shown to be realistic for this type 
of analysis. On the other hand, their assumption of 
constant physical properties is not accurate. The tempera- 
ture difference that typical char forming ablators could 
experience during reentry is of the order of 4000°F to 
5000°F. Over this temperature range change in physical and 
thermodynamic properties is significant. 

Fundamental Equations : The coordinates used by Kondo 

et. al. (8) are shown in Figure 2-7. The y axis is 
coincident with the reverse direction of reaction which 
extends from zero to The origin is taken at the initial 

position of the solid surface at y=0. The equation of 
continuity for the solid material was shown to be: 

"E/RT 

3 (p_ a) = - p a A (e -a) (2-23) 

Tt 5 3 
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where 


-E/RTi 

a = e [-00 (2-24) 

The equation of continuity for the pyrolysis gases is: 

-E/RT 

3 (p (1-a) ) +3 {p u (1-a) } = p a A (e -a) 
Tt y Ty y y s 


(2-25) 


where "a" is an effective cross sectional area of the solid 

I ■ 

material in the pyrolysis zone. It is defined as "the ratio 
between the density of char at the surface and that of the 
virgin material/ since the density of char crushed into 
powders equals to that of virgin material according to the 
experimental results. 'a' changes to 1 on the boundary of 
the virgin material". (8). 

The equation of motion is: 

p g + P g u g 3u g «= - 3P (2-26) 

3 y STy 

t ’ 1 

and the energy equation is: 



+ k g (1-a) 


-E/RT 

-p s A Q (e -a) + 


H] • h [ p ^ Cp ' 9 <1 - > 1" 

5T C P/ S T “j + ft |^g c P-g (1 “ a> 



(2-27) 

The first bracketed term on the left hand side expresses 
the increase in thermal energy by heat conduction. The 
second term on the left hand side represents the thermal 
energy absorbed by the pyrolysis gases due to convection. 
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The right hand side of the equation represents the energy 
change by pyrolysis reaction (first term) and by non- 
steady effect (last two terms) . 

Finally: 


P = 1 p RT (2-28) 

f V g 

vhich is the ideal gas equation of state. Equations (2-23) 
and (2-25) to (2-28) are the fundamental relationships used 
to solve for the system of five unknowns. P, Pg, T, a and 

u g * 

Coordinate Transformation : Kondo, et:. ad. (8) took 

advantage of their quasi-steady assumption and performed a 
coordinate transformation, v, the regression velocity was 
assumed constant and therefore, the coordinate transforma- 
tion is: 


x = y + v t (2-29) 

which means that the coordinate system moves with the 
regression velocity of the pyrolysis zone. 

Performing this transformation the equation of conti- 
nuity for the solid becomes: 

(a v) = - a A (e" E/RT - a) (-30) 

dx 

and that for the gas : 

-E/RT 

~ (Pa (1-a) (v+u ) } =p a A (e -a) 
dx g g s 


(2-31) 
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The equation of motion becomes: 


p g (u g + V) ^ = - dP 

dx dx 

and finally the energy equation is transformed to: 
fk s a+kg (1-a) > ^ + <* s j§ - k,, S| - p„ u„ C, 


(2-32) 
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(l-a) -p s v C p;S a) g +{p g u g C pjg g -p s v C 
■ - p s a A (e-1) 


da}»ji 

P' s dx 


(2-33) 

Combining Equations (2-30) and (2-31) , integrating 
and evaluating the constant of integration, Kondo, et. al. 
(8) obtained the following equation: 


P g u g 88 v( Ps - P g ) 


(2-34) 


Since p s >> p Equation (2-34) becomes: 


Pg u g * Ps v 


(2-35) 


The authors assumed that since u and v are small in 

g 

comparison with the sound velocity, the pressure P can be 
considered constant. Moreover, the authors assumed that 
kg << k s , which is generally the case. By neglecting 
3^ and substituting Equation (2-35) into (2-33) , the 
energy equation was shown to be: 
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(2-36) 
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Equations (2-30) and (2-36) with boundary conditions: 
= T o or (dl) 0 = ? iven at x = 0 

(2-37) 

“ a o p char at x = 0 
^s 


and 


T = T_ oo and a = 1 at x = - « 


(2-38) 


There are four boundary conditions for the first and 
second order differential Equations (2-30) and (2-36) invol- 
ving two unknowns T(x) and a (x) . Therefore, the regression 
velocity v was determined as an eigen-value of the 
eqiations (8) 


Non-Dimensional Equations : Kondo and co-workers (8) 

found it convenient to non-dimensionalize Equations (2-30) 
and (2-36) by defining the following non-dimensional 
groups : 


Q* = Q/C p#g T.„ 

0 = T/T_ 00 

0 = T /T _ 

0 o' “°° 

E * = E/RT.oo 


X/L = K 


3 = AL/ve- E * /G 


(2-39) 

(2-40) 

(2-41) 

(2-42) 

(2-43) 

(2-44) 

(2-45) 


In addition the author chose to define a characteris- 
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tic area L 2 = k g /(a C pjg p g e E */0 o ). 

Substitution of the above terms in Equation (2-30) 
results in 


da 


E */Q -E*/0 

3 a e ° (e ° -a) 


(2-46) 


and into (2-37) results in a dimensional form of the 
energy equation which is: 


d 2 0 

d£ 2 




1-a 


+ 3 e 


E*/0, 


< e 


-E*/0, 


- a)} fr + 


d© 

cH 


E*/0 -E*/0 E*/0_ -E*/Q 

(y-1) e (e -a) 0 + Q* e (e -a) = 0 

1 (2-47) 

while the new boundary conditions are: 

a = a , © = 0 or « ■ given at 5 = o 

o ° dfo 

and 

a = 1, Q* = 1 at ^ = -oo 

Equations (2-46) and (2-47) were solved numerically 
for the above boundary conditions. The value of the para- 
meters y= 1»0, E*=3 . 0 , a 0 =2 3, and T_ oo =300°K were used in 
the numerical solution and were said to be typical of 
phenolic resin composites (8) . 

Numerical Solution : Equation (2-46) and (2-47) have a 
family of solutions depending on the value of 3 and 
(d©/ d£) Q respectively. However, there is a unique value 
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of 8 that satisfies the condition that at C = -°°, a=l. 

The solid line in Figure 2-8 shows this solution, which 
the authors choose to call a: OK. The dotted line repre- 
sents a solution for which the value of 8 does not satisfy 
the final of a=l. Similarly, there is a unique value of 
(d0/d£) which satisfies the condition that at £= -<», 0=1. 
Figure 2-9 illustrates the correct solution, line ©:0K. 

The dotted line is a solution which does not satisfy the 
final value of 0=1. This type of final value problem is 
usually solved by trial and error. However, Kondo, et. al . 
(3) choose to solve graphically for these values of 8 and 
(d0/d£) Q . This graphical solution is illustrated in 
Figure 2-10. The curve a:0K is the locus of the points 

corresponding to the values (d0/d£) and 8 to get the 

o 

solution a (£) , which satisfies the boundary conditions 
for a. The curve 0:OK indicates the combination of (d0/dy) o 
and 8 which gives the solution 0(5) which satisfies the 
boundary conditions for 0. Therefore, the intersection of 
these two curves determine the appropriate values of 
(d0/d£) o and 8 which satisfy simultaneously both equations. 

In Figure 2-11 it is shown that the thermal decompo- 
sition process takes place over a finite region. Note 
that da/d£ given by Equation (2-46) is the non-dimensional 
expression of the decomposition process. As stated by 
Kondo, et. al. (8) many researchers have assumed that the 
pyrolysis reactions occur at the surface of the "heat 
conducting region" (which is the virgin material) and have 
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neglected the depth of the pyrolysis reaction as is shown in 
Figure 2-11. This, the authors claimed, is inadequate. 

The effect of the change in the non-dimensional 
temperature, 0 , on the temperature profile 0, is illus- 
trated in Figure 2-12. 

In summary, the ablation analysis of Kondo and co- 
workers (8) is for a one-dimensional ablation model. The 
most significant contribution of their analysis is to have 
indicated the importance of considering decomposition in 
depth to take place in a final region of space. However, 
their assumption of constant physical properties is 
considered to be a serious limitation in the analysis. In 
addition, the assumption that phenolic resin decomposition 
is described by a single Arrhenius expression is inadequate. 
As has been indicated before, several pseudo-order kinetic 
expressions of the Arrhenius type are required to 
adequately describe decomposition of phenolic nylon resin 
composites (1, 2). 

In the following sections typical transient response 

f 

analyses along with a summary of related work of flow in 
the char zone are discussed. This summary is given to 
support our mathematical analysis and to justify some of 
the assumptions made. We refer to the work of April (19) 
for a more detailed discussion of the subject of transpira- 
tion cooling, carbon deposition and related areas. 
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Transient Analyses to Describe the Response of Thermal 
Protection Systems " ~ 

The need for the accurate description of the thermal 
response of ablative composites is required to insure a 
safe design for thermal protection of reentering space 
vehicles. The necessity is also dictated by the desirabi- 
lity of reducing expensive ground testing to a minimum. 
Although these facilities have improved over the years, 
all the complex interactions of the shock heated gas and 
the ablative heat shield are not easy to simulate in these 
ground test facilities (55, 93-99). Therefore, accurate 
and realistic mathematical analyses are needed to 
effectively predict the thermal response of ablative 
composites under a wide variety of external flow field 
conditions. 

Almost all major aerospace companies and governmental 
agencies interested in ablative heating have analyses of 
the transient response of thermal protective systems. The 
literature on this subject is abundant. Lapple and 
co-workers (106) have compiled over 525 references on the 
area of ablation of related topics. More recently Penner 
and Olfe (87) in their chapter on ablation, have more than 
thirty references on major aspects of the ablation process, 
for both charring and non-charring materials. To try to 
cover just a small percentage of these references would be 
a difficult and time consuming task without contributing 
much insight into the analysis proposed in this research. 
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Moreover, April (19) has written a detailed review on 
some of the most important work in transpiration cooling 
and has an extensive list of references on the more 
important aspects of the ablative process. In addition to 
the three previously mentioned references, there are some 
important review articles on the topic of ablation which 
have extensive bibliographies on the subject (100-106) . 

For these reasons, it would be considered repetitious to 
go over this material which has been thoroughly covered 
in other places. Rather it is best to devote a brief 
discussion to what we consider to be four important and 
typical transient ablation analyses. These are those of: 
Swann, et. al . (58) , Kendall, et. al. (51, 53) and Kratch, 
et. al. (83) , and Clark (111) . In addition we will 
summarize in a later section some. of the research done on 
flow of pyrolysis gases through porous chars and related 
areas. 

The ablative analysis of Swann, et. al . (58) is a 
one-dimensional analysis of the transient response of 
plastic composites. The thermal response of materials was 
described with as many as three different layers, and the 
first two could have moving boundaries. The coupling of 
the flow field to the ablative surface was done by means of 
an energy balance. This energy balance was performed at 
the char surface, where the convective heating rate was 
computed using either a linear or quadratic approximation 
to the blocking effectiveness for a laminar boundary layer. 



56 


The cold-wall convective heating rate and the radiant 
heating rate incident on the surface were specified func- 
tions of time. These values were used as inputs in their 
numerical solution. The surface removal mechanisms were 
considered to be vaporization at the sublimation tempera- 
ture# and diffusion and/or reaction controlled oxidation 
of the carbon at the surface. 

The energy equation for the one-dimensional# non- 
steady flow of pyrolysis gases in the char zone was given 
as: 
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The term in brackets has been labeled the reacting gas 

n 


heat capacity. 


In this bracketed term# 


E H . R . accounts 

j-i 3 3 


for the heat effects of chemical reactions. When the flow 

is frozen# that is# no chemical reactions take place# the 
n 

term E H. R, equals zero. However# when there are 

j-1 3 3 

chemical reactions# this term is non-zero. This type of 
analysis is used to predict the transient one-dimensional 
thermal performance of charring ablator heat protection 
system when exposed to a hyperthermal environment. 

The analysis of Kendall, et. al. (51-53) is a transient 
one-dimensional analysis of the coupled, laminar# chemically 
reacting# boundary layer to the ablative surface. 
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The boundary layer solution was related to the shock layer 
by specifying edge boundary conditions. Similarly, the 
boundary layer solution was related to the transient res- 
ponse of the ablative composite by surface conditions. 

Four options were available to couple the laminar, 
chemically-reacting, boundary layer to the char surface. 
These included specifying 1) wall enthalpy, 2) mass flux of 
each gas species as they entered the char zone, 3) equili- 
brium mass flux at the surface or 4 ) coupled mass and energy 
balance at the wall as provided by a transient charring 
conduction solution. 

The one-dimensional in-depth analysis considered that 
the decomposition of the virgin material could be described 
by means of an equilibrium analysis; that is, the 
decomposition reactions were in chemical equilibriu. The 
assumption of chemical equilibrium between the pyrolysis 
gases and the char, results in an overprediction of carbon 
deposition. To avoid this problem Kendall and co-workers 
(51-53) assumed that no carbon deposition occurred, and 
that only gas phase reactions took place. A modified form 
of Darcy's Law was used to calculate the pressure drop 
across the char layer. 

The analysis of Kratsch, et. al. (83) was a one- 
dimensional, transient analysis of the coupled mass and 
energy balance in char forming ablative composites. 
Depolymerization of the plastic ablative composite was 
modeled by an Arrhenius type kinetic expression obtained 
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from thermogravimetric data. The receding surface boundary 
conditions and the convective and radiative heating were 
specified input functions to the numerical calculations. 

In addition, chemical errosion of the surface was specified. 
The in-depth analysis of the char layer considered the 
pyrolysis gas products to be in thermo-chemical equilibrium. 
However, the assumption of chemical equilibrium between 
the pyrolysis gas and the char, as previously mentioned, 
results in an over prediction of carbon deposition. To 
avoid this problem only gas phase reactions are considered 
and the amount of carbon was empirically adjusted to that 
found by experiments. 

This work was one of the first attempts to describe 
a rather complex system with an analysis that was not res- 
tricted to simplifications of frozen flow, constant 
physical properties, or omitted heat absorption terms. 

Another work that treats the complex physico-chemical 
phenomena of a charring ablator is that of Clark (111) . 
Clark's work parallels that of the already mentioned work 
of Swann (58) , but Clark expands the mathematical treat- 
ment by taking into account, in more detail, various 
thermal, chemical and mass transfer processes present in 
ablation. Clark's work is distinctive from that of Swann 
mainly in three respects. One is a detailed chemical 
kinetic treatment of the pyrolysis gases in the char. The 
other is a method of describing mass deposition in the 
char. Finally, and most important, is the complexity of 
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taking into account thermal non-equilibrium effects 
between the pyrolysis gases and the char. One interesting 
observation is that Clark uses Stroud's (5, 6) reaction 
plane approximation to take into account the thermal de- 
composition process of plastic composites. His stated 
reason is that such an assumption was necessary to provide 
a second boundary condition for the gas-momentum conserva- 
tion equation. What Clark surprisingly fails to mention 
is that the reaction plane approximation to describe the 
decomposition of the virgin plastic is valid only when it 
can be described by one reactable species. This point is 
stressed because Clark's numerical examples concern 
phenolic-nylon which we know can be described by at least 
seven degradation reactions. We know from the work of 
Stroud that the higher the surface recession velocity, 
the better is the reaction plane approximation. However, 
as we mentioned before, Stroud did not quantify this 
velocity, neither did he speculate as to how good his 
reaction plane approximation would be for those materials 
such as phenolic-nylon, where more than one reactable 
species is present. Clark is not clear either in this 
respect. Nowhere in the text of his work was he explicit 
about the limitations of the reaction plane approximation, 
although we have to assume that he was aware of them. 

In summary, Clark's analysis has the following 


salient features: 
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1) One-dimensional transient. 

2) Multilayer system (char layer, uncharred layer 
and insulation layer. 

3) Intereaction of the external surface with the 
boundary layer by taking into account surface 
oxidation, sublimation and mechanical erosion. 

4) Homogenous and heterogenous chemical reactions 
within the char layer which account for mass 
deposition. 

5) Thermal non-equilibrium between the char layer 
and the pyrolysis gases. 

6) Reaction plane approximation to account for 
decomposition of uncharred layer. 

Clark' s work was very thorough and extensive. He 
studied in detail the various thermal, chemical and mass 
transfer phenomena in the ablative process. His analysis 
predicts that the overall performance for a low density 
phenolic-nylon is 16 percent greater them the performance 
calculated by the emalysis of Swemn (58) . Clark attributes 
this difference in predicted performance to the considera- 
tion of char layer deposition in his analysis. 

Although Clark does not state the actual CPU time 
required to solve his model, he does state that the compu- 
ter solutions are time consuming. The objective of his 
program was to develop the capability of analyzing 
ablation systems will all the complicating factors already 
mentioned to provide guidance in selecting the most 
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important effects. In this way, he provided means of 
calibrating less complex analyses to account for those 
effects which are found to be significant. 

The next section presents a summary of previous 
research of flow in the char zone. 

Summary of Previous Research of Flow in the Char Zone 

One of the most recent studies of flow in the char 
zone was done by April (19) . His mathematical analysis 
consisted in a non-equilibrium, one-dimensional, steady 
flow model which accurately predicted the energy transfer 
in the char zone of a nylon-phenolic resin composite, for 
front surface temperatures of up to 3000°F. The important 
chemical reactions and kinetic data for a temperature 
range of 500°F to 3000°F, with experimental simulation to 
2300°F, were determined and incorporated into the mathema- 
tical analysis. His analysis, in conjunction with 
experimental results obtained in a Char Zone Thermal 
Environment Simulator were used to show the shortcomings 
of the limiting cases of frozen and equilibrium flow 
analysis in predicting the true behavior within the char 
layer. Comparisons of the experimental data for low 
density phenolic-nylon chars were made with the results 
obtained using graphite as a simulated char. The non- 
equilibrium flow analysis was used to accurately predict 
the energy transport in the graphite medium using the same 
important reactions and kinetic data developed for flow 
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through chars. April (19) conducted studies to determine 
carbon deposition and decomposition product distribution 
for methane and phenol using carbon-14 tracers. Carbon 
deposition measurements within the char layers were used 
to locate the temperature where chemical reactions among 
the pyrolysis product became significant. In addition, 
oxidation of nylon-phenolic resin chars were studied 
to determine the rate of oxidation of the char with dis- 
tance from the front surface. It was found that oxidation 
was taking place at all depths within the char; that is, 
reaction rate limited. This, however, is not surprising 
because it confirms results published by Scala (109) . A 
graphical illustration is given in Figure 2-13. April's 

i 

experimental data did not exceed 2300°F which falls outside 
the diffusion 1 controlled regime as shown in Figure 2-13, 
and this is consistent with the results presented by 
Scala (109) . 

A considerable amount of information pertaining to 
the problem associated with the formulation of an accurate 
ablative analysis was discussed by April (19) or Scala 
(109). Others that have contributed to the understanding 
of the flow of gases through porous matrices have been 
Koh and del Casal (11, 12), Clark (16), and Weger and co- 
workers (17, 18). These studies have been especially 
useful in evaluating the magnitude of the terms in the 
equation of change in the development of more accurate 
and realistic ablative analyses. For this reason the 
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important conclusions and recommendations drawn from these 
authors will be summarized below. 

Thermal Equilibrium Between the Gas and Char ; It is 
important to determine the existence of thermal equilibrium 
since this can greatly reduce the complexity of the equa- 
tions of change (11) . If a small temperature difference 
exists between the gas and the solid matrix^/, it is likely 
that the effect on the energy transfer will be negligible. 
April (19) measured the temperature of the pyrolysis gases 
1/4 inch away from the back surface of the experimental 
char, and found the difference to be approximately 200°F. 

He reasoned, that by the time these pyrolysis gases would 
arrive at the back surface of the char the thermal gradient 
should be further reduced. Clark (16) reported differences 
of 200° R to 800°R at the midpoint of 0.021-0.33 feet thick 

graphite and carbon matrices over a wide range of mass flux 

2 

values, 0.018 to 0.07 lb/ft -sec. These large differences 
were the result of the initial gradient (2000°R-3000°R) 
between the gas and solid at the matrix back surface. This 
abnormally large gradient was produced by the resistance 
heating apparatus used to simulate high temperature re- 
entry ( 4000°R) . This however, is not representative of 
the conditions encountered on thei back surface of the char 
during reentry. Instead, the pyrolysis gas and char back 
surface temperatures are approximately equal to the plastic 
decomposition temperature ( 1500°R) . Since no large 



65 


initial gradient between the gas and solid phases exists 
this assumption should be valid. 

Variable Physical Properties : The assumption of 

constant physical properties is only valid over a relative- 
ly small temperature range. In ablative cooling where the 
temperature can exceed 3000°F changes in physical proper- 
ties must be expected. j 

Koh and del Casal (11) noted significant differences 
in the results obtained between the constant and variable 
fluid physical properties model. Their findings were for 
solid matrices for porosities greater than 0.5. This was 
apparently attributed to the increasing importance of gas 
convective heat transfer effects in the high porous 
matrix. The only modes of energy transfer accounted for 
were gas convection and solid conduction. In ablative 
cooling using nylon-phenolic resin chars with porosities 
between 0.7 and 0.8, the gas convective term will be 
significant and hence, properties must be considered 
variables over the large temperature range. 

With regard to the solid, or char properties, Weger, 
et . al . (17, 18) measured the change of char porosity 
and permiability when carbon deposition, or depletion 
occurred. They found these properties to have measurable 
changes. These results, however, were obtained for 
materials with porosities of 0.21 to 0.35. The effects of 
these changes on the pressure drop and energy transfer for 
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high porosity (0. 7-0.8) materials were not determined. 

Until more pertinent data is obtained for high porous 
materials, the permiability and porosity have to be taken 
as constant in the analysis. This is a good approximation 
in any porous material where excessive carbon deposition 
or depletion does not occur. 

Modified Form of Darcy's Law ; This equation is used 

to describe the momentum transfer for flow through porous 

media. Darcy's empirical equation relates the fluid 

velocity to the pressure drop within the porous media. 

Weger, et. al. (17, 18) found that a modified form of 

Darcy's Law gave a better prediction of the pressure drop 

across porous materials. This modified form of the law had 

an additional term which accounted for the inertial effects 

These effects become significant for mass fluxes exceeding 
2 

(0.01 lb/ft -sec). For ablative cooling mass fluxes can 

2 

be as high as 0.05 lb/ft -sec. Hence, the modified form 
of Darcy's Law should be used. 

Various Modes of Energy Absorption ; In all the work 
reviewed in this chapter which dealt with the temperature 
distribution in porous media, gas convection and solid 
conduction were considered to be major modes of energy 
transfer within the char. Recently, the Southern Research 
Institute has investigated the effect of internal radiation 
as a mechanism of heat transfer (37) . While admitting that 
this mechanism can be important it was difficult, 
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experimentally, to separate this effect from others. 
Therefore, they incorporated the radiation effect in the 
thermal conductivity they reported. 


Gas conduction in most cases was ignored. Koh and del 

Casal (12) presented results which clearly defined when 

gas conduction could be neglected. This occurred for 

either high mass flow rates, or small temperature 

differences across the porous material. At low flow rates 
2 

( 0.01 lb/ft -sec), or large temperature differences 
( 3000°F) , energy absorption by gas conduction becomes 
important. Since -these conditions are present in ablative 


cooling, gas conduction must be taken into account in any 
realistic analysis. The effects of chemical reactions on 
the energy absorption and the physical properties of the 
material are discussed below. 

Effects of Chemical Reactions : The various effects 

that chemical reactions have on flow through porous media 
have been discussed by a number of authors (13, 16, 17, 18) 
Koh and del Casal (13) and Clark (16) considered chemical 
reactions as important modes of absorption, while Weger, 
et . al . (17, 18) used them to explain the changes in the 

physical properties of the porous specimen due to carbon 
deposition or depletion. In any realistic analysis 
involving high temperature flows, chemical reactions must 
be considered. Clark (16) found in his studies of methane 
flow through porous carbon and graphite that there were 
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three distinct chemically reacting zones. At low tempera- 
tures, a frozen flow region where no significant chemical 
reactions occurred was observed. As the temperature in- 
creased an intermediate region where chemical reactions 
were kinetically controlled was noted. Finally, the 
third distinct region was that in which reactions were in 
chemical equilibrium. It is likely that for more complex 
chemical systems, such as the pyrolysis products produced 
from the degradation of plastic composites, the gases 
would undergo a similar transition from frozen through 
non- equilibrium to equilibrium. Therefore, a non-equili- 
brium flow analysis would be required to accurately 
predict the energy absorption within the char layer. Weger 
et . al. ( 18 ) also observed a non-equilibrium region where 
the reactions were kinetically controlled. 

We have established the importance of a non-equilibrium 
flow analysis as a pre-requisite to an accurate description 
of the energy transfer for ablative composites. Two 
important problems, however, remain to be discussed. One 
is the accurate description of the composition of the 
pyrolysis gases entering the char zone. The other is 
finding the important chemical reactions with thermodyna- 
mically consistent kinetic data, that occur in the tempera- 
ture range of interest. These will be discussed in the two 
remaining sections of this chapter. 

Pyrolysis Gas Composition : In order to realistically 
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describe the flow of pyrolysis gases through the char zone, 
and consequently, obtain an accurate prediction of the 
energy transfer, the pyrolysis gas composition must be 
accurately known. The current state of the knowledge 
precludes prediction of the products for decomposition 
reactions for even relatively simple polymers. Prediction 
of the pyrolysis products from a mixture o fj polymer, such 
as those contemplated for use in heat shields, is 
especially difficult and therefore must be determined 
experimentally. Reactions of the pyrolysis product among 
themselves and with the pyrolyzing resin further com- 
plicates the analysis. 

In all previous research to date, simplified gases 
(helium, methane, etc.) or gas mixtures (methane, carbon 
monoxide, hydrogen, etc.) were used. While these are 
part of the decomposition products, they do not constitute 
the entire product composition expected from the decomposi- 
tion of phenolic-nylon plastic composites. Nelson (1) and 
Sykes (2) addressed themselves to the problem of identi- 
fying the products of pyrolysis. April (19) used these 
studies concurrently with his own studies, and that of 
others to achieve an accurate description of the pyrolysis 
product stream composition entering the back surface of the 
char. A summary of April's (19) results along with 
additional findings of the product distribution is presen- 
ted in Appendix F. 
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Determination of the Important Chemical Reactions and 
Associated Reaction Kinetic Datal To describe the flow 

of pyrolysis gases, and hence, to correctly predict the 

energy absorbed by chemical reactions it is essential that 

the important chemical reactions be established and the 

accurate kinetic data be determined for these reactions. 

This is important, as previously mentioned, because 

i 

chemical reactions do have a profound effect on the energy 
absorption. To achieve this goal, out of the chaotic state 
of the reaction kinetic data in the literature (107) , it 
was necessary to devise orderly and systematic procedures 
to select, and properly screen the reactions and their 
corresponding kinetic data. A continuous search and 
screening technique has been used during this study to up- 
date the reactions included in the analysis. Each reaction 
chosen must be based on the pyrolysis products composition 
initially present, plus the products produced by the further 
cracking of the pyrolysis gases. 

More difficult than establishing the important chemical 
reactions is the task of locating accurate sources of 
kinetic data for the chosen reactions. Many times there is 
a total absence of kinetic data in the literature. Other 
times several sources of data appear, often, conflicting 
(108). In this research the kinetic data has been analyzed 
for consistency with thermodynamic principles. In addition 
the reverse reaction rate constants have been obtained by 
using the chemical equilibrium constant. The technique 
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used to decide which reaction to include in the analysis 
and which data to use is discussed in Chapter V. 

Summary 

In this chapter a detailed presentation of two 
analyses describing the decomposition in-depth of plastic 
composites has been given. These analyses were those of 
Stroud (5, 6) and Kondo and co-workers (8)/ Stroud's 
analysis was mainly concerned with determining whether the 
assumption of a decomposition reaction taking place in a 
plane, could be considered a valid engineering approxima- 
tion. He found that his reaction plane approximation would 
hold for a quasi-steady, one-dimensional flow of gases, 
where one species was considered reacting and where the 
surface recession velocity was high. Kondo and co-workers 
(8) , on the other hand showed the importance of considering 
the decomposition reactions to take place in a finite 
region. As in the case with Stroud they considered only 
one degrading reaction to describe the decomposition of 
the plastic. Both analyses used an Arrhenius type expres- 
sion to describe the kinetics of the decomposition process 
of the plastic. It was established that Arrhenius type 
expressions were independent of the heating rate and thus- 
suitable for these analyses where the heating rates could 
vary . 

The transient analyses of Swann, et. al. (58) , Kendall 
et. al. (51,53), Kratch, et. al. (83) and Clark (111), were 
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briefly discussed. In addition, a brief summary of previous 
research of flow in the char zone was presented. These 
were the works of Koh and del Casal (11-13) , Weger, et. al . 
(17, 18) , and Clark (16) , which showed the importance of 
considering chemical reactions as significant modes of 
energy transfer and as mechanisms which affect the physical 
properties of the flowing gases and the solid matrix. These 
and other considerations have been taken into account in 
the analysis for describing both decomposition in-depth and 
the flow of gases in the char zone, in this research. This 
will be shown in the next chapter. 
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CHAPTER III 


DEVELOPMENT OF THE MATHEMATICAL ANALYSIS FOR 
THE EVALUATION OF THE ENERGY TRANSFER IN 
THE DECOMPOSITION ZONE AND CHAR 
ZONE OF A CHARRING ABLATOR 


Introduction 

I 

The previous mathematical analyses (1, ; 2) to pre- 

dict the energy transfer in the char zone are handicapped 
by the limitation that the char zone depth has to be 
specified a priori . At the point selected as the inter- 
face between the char and pyrolysis zones , it is necessary 
to specify the interface temperature, heat of pyrolysis, 
pyrolysis gas mass flux, and pyrolysis gas composition. 
However, physically there is not a sharp interface between 
the two regions since polymer degradation, pyrolysis gas 
generation and char formation occur over a rather wide 
temperature range. In fact, polymer degradation begins 
at about 200°C, and about 80 percent is pyro lysed at 475°C 
(3) . At this temperature heat absorption by the pyrolysis 
gases is becoming significant as a result of changes in 
sensible enthalpy and chemical reactions. Consequently, 
the combined energy absorption must be predicted from 
considering both the polymer degradation and the gas-char 
interaction simultaneously. 
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The pyrolysis products, formed by the thermal 
degradation of the plastic composite enter the char layer 
and transpire through the front surface of the char, at 
Z = L. A schematic diagram illustrates this in Figure 3-1. 
These gases entering from the decomposition zone experience 
a temperature increase as they flow through the porous char, 
and absorb heat by virtue of their own heat capacity 
and due to mostly endothermic chemical reactions. These 
endothermic reactions are due to the further cracking of 
the pyrolysis gases into lower molecular weight species. 
These species react with each other and with the carbon- 
aceous char layer producing further endothermic reactions. 

The description of the momentum, energy, and 
mass transfer equations for the combined decomposition- 
char zone analysis is obtained by simplification of the 
general equations of change (continuity, momentum and 
energy) to forms applicable to a chemically reacting flow 
(equilibrium and non-equilibrium) , through the decomposition 
zone and char layer of a char forming ablator. Finally, 
typical boundary conditions are specified, followed by a 
discussion of the numerical solution of the equations. 

Statement of the Problem 

To predict the energy transferred in the combined 
decomposition-char zone the equations of change are written 
to apply to any point in the flow field. Using a quasi- 
steady analysis, or steady state approximation, the point 
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Figure 3-1. A Sketch of the Virgin Plastic, 
Decomposition and Char Zone. 
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of view is taken of an observer moving with the negative 
surface recession velocity, v. This is illustrated in 
Figure 3-1 which shows schematically the mass fluxes of 

virgin plastic, pyrolysis gases and solid., The point at 

/ 

which thermal degradation of the plastic is initiated is 
taken to be that of Z=0. 

The material balance relating the virgin plastic 
flow with the flow of pyrolysis gases and degrading solid 
in the combined decomposition-char zone as written by 
Stroud (4) is: 

p Q v = W + pv (3-1) 

where W is the total mass flux, v is the surface recession 
velocity, and p Q and p are the virgin and degrading solid 
densities respectively. This assumes that the bulk volume 
of the char is the same as that of the virgin plastic, 
which is generally the case. 

Referring to Equation (3-1) for a known surface 
recession velocity and densities of the virgin and degrading 
composite, the gas mass flux can be computed. In addition, 
the composition of the gases generated by the decomposition 
of the virgin plastic composites must be known to be able 
to accurately predict the energy absorbed due to chemical 
reactions in the char zone. As illustrated in Figure 3-1, 



these pyrolysis products enter the char zone and exit 
through the front surface of the char, at Z=L. Changes in 
the mass flux of the various species within the char 
occur as a result of chemical reactions at finite reaction 
rates, Rj . / 

The particular restrictions and assumptions made 
in the formation of the combined zone analysis are pre- 
sented and justified subsequently. The simplifications 
of the general equations of change resulting from these 
restrictions follow. Finally, the solutions of the re- 
sulting equations for frozen (no chemical reactions) , 
equilibrium (calculated from thermodynamics) , and non- 
equilibrium (calculated from reaction kinetics) flow are 
given. 

Restrictions to the General Equations of Change for Flow 
in the Combined Decomposition-Char Zone 

Several restrictions and assumptions are applied 
to simplify the general equations of change to reduce 
their complexities. In the subsequent discussion, these 
assumptions and restrictions are analyzed and justified. 

Quasi-Steady and One Dimensional Flow of Pyrolysis 
Gases : The assumption of steady flow of the pyrolysis 

gases is based on the fact that after a brief initial 
period of time the thickness of the char has been shown 
to remain constant. The data of Peters and Waddin (5) 
in Figure 3-2 for a 50 : 50 weight ratio of nylon-phenolic 










91 


ablative composite, formed in a subsonic electric arc 
jet shows this graphically. This implies that the 
rate of formation of char, formed from the degradation 
products of the virgin material equals the rate of 
consumption of the char by air oxidation and/or mechanical 
erosion by the flow in the boundary layer. In addition, 
to this, the quasi-steady analysis is justifiable because 
the residence time of a fluid particle of pyrolysis gas 
is small (less than 0.01 seconds) when compared to the 
rate of change of the char surface. 

The one dimensional assumption is easily 
justified since the radii of curvature of reentering 
spacecrafts (e.g., Mercury, Gemini, and Apollo) are large 
(about 5 to 8 feet) in comparison with the char thickness 
(<^1/4 of an inch) and as such, the flow is one dimension- 
al and normal to the front surface. 

Pyrolysis Products Behave as a Perfect Gas Mixture: 
The assumption that the pyrolysis products behave as a 
perfect gas mixture implies that the equations of state 
for ideal gases are applicable. This is a realistic 
assumption considering the high temperatures ( 1 500 - 6000°F) , 
and low pressures (''VI atm) encountered in planetary 
reentry by blunt body vehicles. However, it also implies 
that the heat capacity is that of a real gas and hence, a 
function of temperature. In Table 3-1 the values of the 
compressibility factor, Z, are presented for some of the 
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most important constituents of the pyrolysis gases. The 
temperature at which Z has been calculated is 450°K. This 
temperature is the lowest at which the pyrolysis gases 
are formed. It is also the temperature at which degradation 
of the plastic is initiated. T c and P c ar ; e the critical 
temperature and pressure, and T r and P r are the reduced 
temperature and pressure for the respective gases. As is 
shown in Table 3-1 all of the compressibility factors 
are greater than 0.995; which justifies the assumption of 
ideality. 

Virgin Material, Char and Gas Physical Properties ; 
Because of the existence of a large temperature gradient 
from the front surface of the char layer to the virgin 
material, physical property variations with temperature 
must be considered. Variations with temperature of the 
virgin material's thermal conductivity is accounted for. 
Thermal conductivity variations with temperature have 
been reported by Wilson (5) , Engelke (6) and Nagler (7) . 

For the char, the porosity change with temperature is 
not accounted for due to a lack of experimental data on 
the subject. The value of the porosity is thus kept con- 
stant throughout the analysis. The thermal conductivity 
variations with temperature of the virgin material, the 
pyrolysis gases and the char are taken into consideration 
in the analysis. The data used is that reported by the 
Southern Research Institute (6) and the Jet Propulsion 
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Laboratory (7) of the California Institute of Technology. 
These data were fitted empirically and are presented in 
Appendix C. 


Thermal Equilibrium Between the Pyrolysis Gas 

I 

Products and the Char : The temperature gradient between 

the pyrolysis gas products and the porous carbon matrix 
is assumed to be small. This is supported by data pre- 
sented by Koh and del Casal (8) in which they showed a 
maximum temperature difference of 300°F for the flow of 
air and helium through packed beds of spheres. The maxi- 
mum temperature of the matrix was 2700°F. 

Clark (9) also addressed himself to the problem 
of quantitatively determining the existence of a thermal 
gradient. He reported differences of 200 to 800°R between 
the gas and the porous matrix over a wide range of mass 
fluxes, 0.018 to 0.07 lbs/ft^ sec. In this particular 
investigation thermal non-equilibrium between the gas and 
solid is primarily caused by the resistance heating 
method for achieving temperatures between 3000°R and 4000°R. 
In ablative heat protection applications the gas flowing 
from the pyrolysis zone into the char layer is at the 
local char temperature, and the abnormally large tempera- 
ture difference characteristic of the resistence heating 
apparatus is non-existent. Therefore, the assumption of 
thermal equilibrium between the gas and char in ablative 



cooling should be a very good approximation to the real 
behavior . 


Momentum Transfer in the Char Zone : The modified 

form of Darcy's Law was used to determine the pressuee 

/ 

drop across the char layer. This was based on the work 
of Weger et al . (10) . It was shown in this work that the 
inclusion of an inertial term gave a more accurate 
prediction of the experimental pressure drop for values 
of the mass flux of the order of 0.05 lb/ft^-sec. 

PV Work and Viscous Dissipation : The pressure 

drop across one quarter inch thick, low density nylon- 
phenolic resin char was experiementally determined by 
April (2) which. showed it to be approximately 15 lbs/ft 6 
for a pyrolysis gas mass flux of 0.05 lb/ft -sec and a 
front surface termpature of 2000°F. In addition, he 
computed the PV work contribution to the energy transport 
equation and showed it to be 1.2 BTU/ft 3 sec. This was 
compared to the convective energy term, evaluated at the 
back surface where the temperature was small, and was 
shown to be 1000 BTU/ft^-sec for a gradient of 
approximately 40,000°F/ft and an average heat capacity 
of 0.5 BTU/lb-°F. This clearly demonstrated that the PV 
work term in the energy equation was negligible. Since 
the velocity 5 ft/sec) and the viscocity ( /N ^ 0.05 cp) 
of the gas mixture is small, the energy generated by 



viscous dissipation can be neglected also and the term 
omitted from the energy equation. This allows for the 
uncoupling of the momentum and energy equation. Thus 
simplifying considerably the numerical solution of the 
problem. j 

Diffusional Transport : Energy, or mass transport 

by diffusion, is negligibly small in comparison with the 
bulk fluid transport. The average residence time of a 
gas particle in a one-quarter inch thick char layer is 
0.01 seconds for a mass flux of 0.05 lb/ft 2 -sec. 

Derivation of the Equations of Change for Flow in 
the Combined Decomposition-Char Zones 

The application of the above restriction to the 
general equations for flow of the pyrolysis gases in the 
combined decomposition-char zone is now discussed. 

Specie Continuity Equation : Referring to Figure 

3-1 the specie's continuity equation for the ith component 
of a gas mixture for flow through a porous medium is (11): 

Doi _ _ _ _ 

- - Pi (V-u) - (VJi) + Ri (3-2) 

where p-^ is the concentration, Ji, the mass flux by 
diffusion, R j , the rate of generation of chemical specie i 
and u is the velocity of the pyrolysis products within 
the pores. 
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For a steady, one dimensional flow of gases, 
neglecting mass transport by diffusion, Equation (3-2) 
reduces to: 

3 ! (Piu) = Ri i (3-3) 

Assuming that the pyrolysis gases do not lose or 
gain any mass by reaction. Equation (3-3) becomes, summing 
over all the gas species: 

l af ( p i u) = ° (3-4) 

i= 1 

However, if the pyrolysis gases lose weight due 
to carbon deposition. Equation (3-4) becomes: 

I at ( p i = - R c (3-s) 

i = l az 

where R c is the amount of carbon deposition. If we de- 
fine Wp as the mass flux of pyrolysis gases based on the 
cross-sectional area of voids in the char (units of 
lbs/ft voids'sec) we can say that: 

dWp n i 

“HI- = l ml d? ( pi u) - - Rc 

Equation (3-6) says that any change of the mass 
flux of pyrolysis gases is due to the loss of carbon. 
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Momentum Equation : The momentum equation for flow 

through porous media was formulated by H. P. G. Darcy in 
1856 (8). Darcy observed during experiments with a one 
dimensional packed bed that gas velocity at any point in 
the bed was directly proportional and in the same direction 
to the pressure gradient at that point. In vector notation, 
including the effect of body forces when considering a 
vertical flow direction, Darcy's Law is: 

u = - • (VP - Pg) (3-7) 

«» 

Applying this equation to a one-dimensional, horizontal 
flow through a porous char layer and solving for the 
pressure gradient gives: 

- 4 ? = Y (UE) (3 ' 8) 

This equation is valid at low gas flow velocities 
within the porous medium. However, at high gas velocities 
it is necessary to add a term to account for the inertial 
effects. This additional term leads to a modified form 
of Darcy' s Law: 


dP 

- ar = 


~ (ue) + Bp (ue) 2 


(3-9) 



Multiplying both sides of Equation (3-9) by the 
gas density, p, followed by substitution of the ideal gas 
equation of state (p = PM W /RT) on the left hand side of 
the Equation (3-9), results in Equation (3-10): 

/ 


V dp 

- IT * ar = P Cm/30 


(ue) + Bp 2 


(ue) 2 


(3-10) 


If we define W as the total mass flux of pyrolysis 
gases based on the total area, we have that: 

W = eWp = epu (3-11) 

therefore substitution of Equation (3-11) into Equation 
(3-10) results after rearrangement: 

- PdP = jp [ (M/ Y) w + 3 (W) 2 ] (3-12) 

w 

Integration of Equation (3-12) between the front 
surface pressure (P = at Z = L) and any point within 
the char layer, (P at Z) , results in an integral equation 
for the pressure distribution over the char. 


L L 

P = (P L 2 + 2R [/(m/y) IV (T/M w )dz + / 3 (T/M w ) 
z z 


(W) 2 ds3 h 

(3-13) 


In this equation all parameters that vary with 
temperature (hence, char distance) are left under the 
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integral signs. These variations are calculated by 
polynomials in temperature and from the solution of the 
energy equation. 


Energy Equation . To formulate the^ equation that 
describes the energy transfer in the combined decomposi- 
tion-char zone, the energy equations for the gas and 
for the solid are written separately. They are later 
combined considering that the gas and solid are at the 
same temperature at any cross section in the flow; that 
is in thermal equilibrium. 

The general form of the energy equation for a gas 
mixture containing n species is (11) : 

r^'T __ ___ n _ 

PC Sr = - (V-q) - (T:VU) + l (Ji- gi ) + 

P uz i+1 

( t-itt) bt * | =1 «i icv-Ji) - Ri) (3-i4) 


For a one dimensional steady flow of gases in the 
combined zones, neglecting viscous dissipation, work 
against gravity and diffusional effects, compared to the 
heat transferred by ''.onduction , convection and chemical 
reactions, Equation (3-14) becomes: 




5 In V 
d In T 


udP 


n 


I Hi R. 
i = l 


(3-15) 



For an ideal gas — in T * s one CH) • Further- 
more, the work by pressure forces across a high porous 
char can be neglected as has been shown in the previous 
section and Equation (3-15) simplifies to: < 

pC P U 37 = - 4 (qz) -l »i Ri (3-16) 

dz dz i _2 

where q z is, from Fourier's Law of heat conduction: 


dT 

<*z = ’ k g 37 


(3-17) 


The energy equation for a reacting gas flow 
through a porous media follows directly from multiplying 
Equation (3-16) by the porosity e and by substituting 
the definition of Equation (3-11) and Equation (3-17) 
into Equation (3-16) and is: 


e * F * C P 37 = e HI t k g 37 * ' E | =1 Hi Ri (3 ‘ 18) 

where e • jr • Cp ^ is the convective term 


(kg ^ ) is the conduction term 
n 

e I Hi accounts for the thermal effects of the 
i = l 

gaseous chemical reactions. 
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Similarly, the energy equation for the solid in 
the combined decomposition-char zone simplified for the 
same restrictions to the following form: 


/ 


(i-eKp^ v 27 = ( 1 -e) 27 ( k s §7^ + C 1 ’ 6 ) l_ n+1 Hi 


dT 


Ri 


q(T) 


where 


(3-19) 


( 1 -e) C^id V ^7 i s the convective term 
( 1 -e) 27 ( k s 37 ) is tke conductive term 
1 

( 1 -e) £ Hi Ri accounts for the thermal effects 
i=n+l 

due to chemical reactions of the 
char. 

q(T) accounts for the thermal effects due to degradation 
of polymer. 


The q(T) function for this research has been 
obtained from data reported by Sykes (12) for important 
phenolic-nylon and silicone elastomers heat shield 
components. Typical curves reported by Sykes for 
phenolic micro-balloons, phenolic resin and nylon are 
shown in Figure 3-3. 

The total energy transferred in the combined 
decomposition-char zone is formed by adding Equations 
(3-18) and (3-19) and by using the definition of W of 
Equation (3-11): 



Kj / Kg 


AT = 0 


NYLON 
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TEMPERATURE, ° C 


Figure 3-3. Differential Thermal Analysis 
Thermogram of Nylon, Phenolic 
and Phenolic Microballoons as 
Reported by Sykes and Nelson, 
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It C p Wp * (1-e) C Pj 9 vl § - 3 ! (k e g)- 


1 

l H i R i + q(T) (3-20) 

i = 1 

/ ■ 

where k e represents an effective thermal 
conductivity defined as: 


k e = e + (1-e) k c (3-21) 

1c is the thermal conductivity of the flowing pyrolysis 

o 

gases and k that of the solid matrix. 

1 

In addition the term l Hi Ri is defined as: 

i= 1 


1 

I Ri 

i=l 

The above represents the energy absorbed by the 
chemical reactions on a "total volume" basis. The 
solution of the energy equation, Equation (3-20), gives 
the temperature distribution in the combined decomposition 
char zone. 

In addition to the equations of change just 
developed, one additional equation is considered. This 
is the equation for the net heat absorption within the 
combined zones. 


n 1 

= l Hi Ri + (1-e) l H R (3-22) 
i=l i=n+l i i 



Heat Flux Equation . The heat flux equation is 
used to determine the net heat transfer to the char 


layer and decomposition zone in the charring ablator. 

It is defined as the difference in the heat flux values 

/ 

at 1=0 located in the virgin material, and that of the 
front surface of the char, at Z=L, and is: 


q net q L ‘ q 0 


^e 


dT 

dz 


z = L 


v dT 

^e HT 


(3-23) 

z = 0 


Solving for k e • ^ term of the energy equation. 
Equation (3-20) , substituting it into Equation (3-23) 
and integrating gives the equation needed to evaluate 
the net heat flux within the combined zones; which is: 



T 0 is the temperature at Z=0, Tp is the maximum 
temperature of the decomposition zone, and Tp is the 
temperature at Z=L. 
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I f we define Q as : 


td 


Q - 


J f#k-<lT * Jd-O C ps p s v dl 

To T o / 


(3-25) 


we have: 


T L 


°lNET 


’ J | £ c Pi W P X i dT * ^ J C Ps p s 


t d 



Hi ffi 
dT/dz 


• dT + Q 


(3-26) 


where the first term of Equation (3-26) represents the 
heat absorbed due to the sensible enthalpy change of the 
gases. The second term accounts for the sensible 
enthalpy change of the solids. The subsequent term 
accounts for the heat absorbed by the chemical reactions. 
This term is calculated under three chemically distinct 
flow conditions. One is frozen in which case the term 
is zero. A second is equilibrium, in which the rate of 
heat absorption is calculated by considering the species 
to be in chemical equilibrium, and, the third is 
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non-equilibrium in which this rate is calculated considering 
the species to be reacting at a finite rate. The finite 
rate or non-equilibrium analysis is the most difficult 
to implement for several reasons. One is the problem 
of selecting the important chemical reactions with 
which to model the chemically reacting flow system. 

A second difficulty is the evaluation of the wealth 
of kinetic data for the important chemical reaction, 
and the selection of that data which will most typically 
represent a chosen reaction. 

Of the two reacting flow analyses, the non- 
equilibrium flow is more complex than the chemical 
equilibrium analysis. The latter involves the solution 
of a set of algebraic equations with the energy equation, 
while the former requires the solution of the energy 
equation with coupled, ordinary non-linear differential 
equations which are the species continuity equations. 

There is sometimes associated with the solution of 
ordinary differential equations, with widely separated 
eigen values, a phenomena called stiffness. This 
phenomena is present in chemically reacting flow systems 
where very fast reactions occur. A more thorough 
discussion of this problem is presented in Chapter V, 
and it will be shown how this phenomena can cause 
numerical instabilities and require the integration step 



size to be very small. Finally, Q, the last term in 
Equation (3-26) is the rate of heat absorbed in the 
pyrolysis of decomposition zone, and was defined by 
Equation (3-25) . , 

A summary of the important differential equations 
for describing the energy transferred in the 
decomposition zone and char zone is given in Table 3-2. 

In the following sections the boundary conditions and 
the numerical solution of these equations are discussed 
which will complete the theoretical discussion of the 
combined zone analysis. 

Boundary Conditions for the Solution of the Energy 
and Species Continuity Equations in the Combined - 
Decomposition Zone and Char Layer ~~ 

The energy equation describing the one dimensional 
flow of heat in the decomposition zone and the char 
layer of a charring ablator is a second order, non- 
linear differential equation with variable co- 
efficients. The energy equation is coupled to the species 
continuity equations which are first order differential 
equations. The momentum equation, as previously 
mentioned, was uncoupled from the energy equation 
because PV work was negligible. 

April (2) , who solved the energy equation for 
the char layer only, considered two sets of boundary 
conditions. The first set specified the temperature 
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at the front surface, and, the temperature and 
pyrolysis gas composition at the back surface of the 
char. These conditions made the solution of the energy 

equation a two point boundary value problem which 

/ 

required an iterative procedure. A second set of 
boundary conditions specified the temperature, 
composition and back surface heat flux and made this 
an initial value problem. For both the initial and 
final value problems the thickness of the char was 
considered to be 1/4 inch. The pyrolysis gas mass 
flux, W, in addition, was considered a parameter 
and was varied between 0.01 and 0.05 (lbm/ft 2 - sec). 

In the present analysis the two point boundary 
value problem was not considered. Rather than 
specifying, in addition to the back surface boundary 
conditions, the char thickness and the front surface 
temperature, the thickness of the char was allowed to 
be determined by the front surface temperature. The 
front surface temperature was a parameter in the 
solution and thus eliminated the necessity for an 
iterative solution. Not only does this simplify the 
problem but it eliminates, in addition, some of the 
artificiality associated with fixing the thickness 
of the char. 



In this analysis the back surface temperature, 
the back surface heat flux and the composition of the 
plastic composite were specified. Unfortunately, the 
current state of the art precludes prediction of the 
composition of the pyrolysis product from (a known 
combination of plastic composites. Therefore, it was 
necessary to specify, in addition, the pyrolysis gas 
composition. How this composition was arrived at is 
discussed in Appendix G. The boundary conditions 
selected were therefore: 


T = 

dT I 
dz 


To 



with Tl as a parameter 



x i x i o » 


i = 1 


n 


(3- 


Once the energy equation is solved, the 
momentum equation is solved by specifying the front 
surface pressure. Because the pressure drop is small 
(-15 lbs/ft 2 ) , the assumption of constant pressure in 
the solution of the energy equation can be considered 
valid. 

The computer program implementation for the 
combined char layer and decomposition in depth is 
presented in detail in Appendix A together with a 
complete listing of the program. 
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As mentioned in the Introduction, this research 
deals with the development of a mathematical analysis 
to predict the energy transfer for the combined 
decomposition zone and char layer. In addition, it 
deals with developing the analysis for the char zone 
to 5500°F. Decomposition in depth is now analyzed, 
and this is followed by a brief description of the 
application of the transport equations to frozen, 
equilibrium and non-equilibrium flow in the combined 
zones . 


Analysis of Decomposition in Depth . It was 
previously stated that the previous mathematical analysis 
(2, 12) to predict the energy transfer in the char 
zone is handicapped by the arbitrary separation of the 
decomposition zone and char layer. This analysis was 
useful in bounding the energy absorption in the char 
zone. However, the energy transferred in the two 
zones can be predicted accurately by a combined 
analysis of both the decomposition zone and char zone. 

A schematic diagram was shown in Figure 3-1 for 
the quasi-steady model with the combined zones. The 
point of view was taken of an observer moving with the 
negative surface recession velocity, v. It was shown 
by a simple material balance (Equation (3-1)) that by 
knowing the density of the virgin and degrading materials, 
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the mass flux W of the gases could be calculated. This 
differs from the previous analysis (2, 12) where the 
mass flux of pyrolysis gases, W, was a specified 
parameter at the back surface of the char. The present 
analysis also differs in that the thickness of the char 
is treated as a variable rather than a specified para- 
meter. 

The analysis of decomposition in depth requires 
that the rate of mass loss, or rate of density change 
with temperature be known accurately. Data for the 
densities of virgin and degrading composites as a 
function of temperature have been reported by Sykes and 
Nelson (3) and Madorski (14). The data of Sykes and 
Nelson (3) is particularly useful since it is for 
phenolic-nylon resins; this data was obtained using 
thermogravimetric analyses techniques. In Chapter II, 
the importance of using a kinetic equation of the Arrhenius 
type to correlate the experimental mass loss rate data 
for polymers, was discussed. It was indicated in this 
chapter that the mass loss of a material was 
affected by the heating rate, and that this effect was 
a source of difficulties for earlier researchers 
modeling the decomposition process. Sykes and Nelson 
(3) used a pseudo-order kinetic expression of the 
Arrhenius type to eliminate this influence. The equation 
is : 
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-dpi 

dt 


l ,o 


[(Pi 


ni 

p c,i)/Pi, 0 } A i ex P [ -Ei/RT] 


(3-28) 


Equation (3-28) expresses the rate/ of change of 
density of a polymer with temperature. Pi,o 
initial or virgin density. Pi is the density of the 
material at temp-rature T, p c the residual or char 
density, and Ai and Ei are the well known frequency 
factor and activation energy parameters. 

In the present mathematical analysis it has been 
assumed that when an ablative composite degrades, it 
degrades independently of the other components; that 
is, no interaction is assumed to occur among the 
composites. There is no known experimental data, at 
least to the author, that takes into account these 
interactive effects. Therefore, of necessity, this 
simplifying assumption has been made. This is 
expressed mathematically: 


dp . p dpj 

at '4 , ~dt 
1 = 1 


(3-29) 


and for a quasi-steady flow the time dependent term can 
be modified by: 


dp - v dp 
dt dz 


(3-30) 



Thus, knowing the kinetic parameters for the 
various components in a blend of virgin material, 
specifying the surface recession velocity, and the 
temperature history from the energy equation, the 
variation in density of the virgin material can be 
predicted by the use of Equations (3-28) , (3-29) and 
(3-30). The pyrolysis gas mass flux is then calculated 
by Equation (3-1) . 

Application of the Transport Equations to Frozen , 
Equilibrium, and Non-kqullibrium Flow of Gases 
in the Combined Decomposition Char Zone 

There are two limiting cases currently used to 
simplify the analysis of the flow of pyrolysis gases. 
These are frozen and equilibrium flow. Frozen flow 
considers that there are no chemical reactions among the 
pyrolysis gases and thus, gives a lower bound of the 
energy absorbed by these gases in the combined 
decomposition-char zone. Equilibrium flow, on the 
other hand, assumes that there are chemical reactions 
which are in chemical equilibrium. This case gives 
the upper limit on the energy absorbed by the pyrolysis 
gases . 

The correct answer, however, lies between these 
two limiting cases. A major portion of this research 
has dealt with the development of a third analysis, 
a non-equilibrium flow analysis which predicts more 
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accurately the actual behavior within the decomposition- 
char zone. In this section the equations of continuity, 
momentum, energy and surface heat transfer are 
applied to develop each of the three flow analyses. 

In the two subsequent chapters the solution of the 
particular equations for each analysis is shown and 
compared with each other. 

Frozen Flow . In the frozen flow analysis the 
pyrolysis products entering the char layer are assumed 
not to undergo chemical reactions. Therefore, the 
composition of these gases remain constant throughout. 

The only energy absorbed by the gases is convective 
energy. This energy is absorbed as the gases transpire 
through the char. This analysis gives the lower limit 
on the amount of energy absorbed in the combined 
zones since it does not account for the heat absorbed 
by endothermic chemical reactions. Of all three cases, 
this is the simplest to solve because the heat absorption 
term due to chemical reaction is zero. This is: 

1 

I Hi Ri = 0 (3-31) 

i= 1 

which simplifies the energy and heat flux equations. 
Applying this to the equations of change previously 
developed results in the following simplifications: 
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Continuity Equation: 

W = Wp e = epu = constant (3-11) 

Momentum Equation: / 

r 2 L 

p = ♦ 2R • C J (u/y) W (T/M w ) dz 

z 

L 1/2 

+ J B (T/Mw) W 2 dzl } (3-13) 

Z 

Energy Equation: 

w ‘ dT d /, dT\ 

U r p w p + C 1_e ^ c p P V 1 37 = 37 ( ke 37 ) ' q(T) 

(3-32) 
(3-33) 

The numerical solution of these equations is 


Heat Flux Equation: 
n 

^NET 


* L 

= l i J e W p C p . dT ♦ q (T) 


discussed in a later section. 
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Equilibrium Flow . The equilibrium flow analysis 
assumes that the chemical species are in thermochemical 
equilibrium. This analysis gives the upper limit on the 
amount of heat absorbed in the combined zojes . The 
reason is that reactions occurring both in the 
decomposition and char zone are predominantly endothermic. 
The set of equations which are applicable to this 
analysis is the same as the equations previously developed: 
Continuity (3-11), momentum (3-13), energy (3-20) and 
heat flux (3-24). But in this case: 

1 

l Hi f 0 (3-34) 

i = l 

The species continuity equation, (3-2), can be 
rewritten in terms of the molal flux of the species, 
and is: 


R i - af ( pu) = 3T (* x O ai (3 ' 35) 

1 

Therefore, in order to evaluate the term \ 

i = l 

which accounts for the heat absorbed by the reactions, the 
mass flux W, and the mass fraction x^ of the species must 
be known as a function of temperature. The species 
composition and molal ratio of gases to carbon are a 
function of temperature* pressure and elemental composition 
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of the virgin plastic and can be calculated by one of 
many approaches (see Chapter IV) . In this study free 
energy minimization is used. 

Non-Equilibrium Flow . This analysis is the most 

/ 

complex of the three. It involves the solution of a 
second-order, non-linear differential equation with 
variable coefficient (the energy equation), coupled to 
n (where n is the number of species) ordinary, non- 
linear, differential equations (species continuity equations. 
For the chemical equilibrium case, on the other hand, a 
set of non-linear algebraic equations were solved coupled 
to the energy equation. Numerically, this is a much 
simpler problem. 

For the non-equilibrium analysis the rate of 
reaction is calculated using finite rate chemistry. 

This requires knowledge of the specific reactions taking 
place, and also requires evaluation of the kinetic data 
for particular reactions chosen. The topic of kinetic 
data evaluation is covered in Appendix D. Once the rate 
of reaction, , of each specie is calculated, the heat 
absorbed by the chemical reaction is treated in the same 
manner as was in the chemical equilibrium analysis case. 

The numerical solution of the energy equation, which is 
the subject of the next section, is essentially the 
same for both analyses. 
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Numerical Solution of the Equations of Change 

It was shown earlier in the chapter that to 
formulate the differential equations to describe the energy 

transfer in the combined zones, the continuity and energy 

/ 

equations, for the gas and solids, were written separately. 
Then these equations were combined by considering that 
there were no thermal gradients between them; the 
result was Equation (3-20): 

cLT d dT \ 

[e C p W p + (1-e) Cj> g v] 37 = 37 (k e jj- ) 

1 

- I Hi Ri ♦ q(T) (3-20) 

i= 1 

This equation is an ordinary, second order, 
non-linear (since the terms within the brackets are 
functions of temperature) differential equation with 
variable coefficient. This type of equation requires 
a numerical integration scheme because there are no 
known analytical solutions for this type of non-linear 
equation. A fourth order Runge-Kutta formulae was used 
because of the accuracy and straight forward nature of 
the self starting method. 
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Numerical Solution of the Differential Energy 
Equation 

To solve Equation (3-20) it is necessary to 
transform the equation to a form suitable for numerical 
integration. This is done by expanding Equation (3-20) 
and solving for the second order term. The result is: 


d 2 T 

dz 2 


[e C p W + (1-e) C p p v 


dk e 1 dT 

-arj* hi 


+ I Hi Ri/(dT/dZ) + q/(dT/dZ)]- jZ (3*36) 

i=i 

Equation (3-36) has the form: 

-- t (T. T) S (3-37) 

Where A represents the terms within the bracket of 
Equation (3-36) divided by the effective thermal 
conductivity, k e . Notice that A is a function of both 
temperature T and the gradient T (= ijX) . 

A commonly employed procedure for the Runge- 
Kutta is to convert the second order Equation (3-37) 
into two first order equations to be solved simultaneously 
(16) . The procedure is as follows: 
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Let 

T = g (3-38) 

Substituting the above in Equation (3-37) gives: 

/ 

« 

= f(T) (3-39) 


Solving Equations 3-38 and 3-39 simultaneously 
gives the wanted solution. A summary of these equations 
is given in Table 3-3. 

The general formulae, as they apply to the 
solution of the differential energy equation is given 
in Table 3-4. The truncation error of this technique 
is of the order O(h^) where h is the step size (16). 

Note that the term f(T,T) of Equation (3-36) 
has an implicit dependence on the mole flux of each 
specie i. This implicit dependence is in the term 
1 

J Hi Ri . Therefore, the species continuity equation 
i= 1 

must be solved simultaneously with the energy equation 
when reactions occur in the flow field. 

Check Procedure. The Runge-Kutta computer 


implemented solution was checked by numerically solving 
Equation (3-37) with f, constant, and comparing it to 
the analytically exact solution of the equation. 
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TABLE 3-3- Conversion of a Second Order Differential 

Equation to the Equivalent set of Two First 
Order Differential Equations. / 



= f(T,T) dT 
dZ 


(3-37) 


which decomposes into two first order equations 
shown below: 


31 = T = G(T) (3-38) 

31 = f(T)= F(T,T) 


(3-39) 
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TABLE 3-4. General Formulation for the Evaluation of 
the Runge-Kutta Parameters for the 
Simultaneous Solution of Two First Order 
Differential Equations (16) . , 


t N + T = T n ♦ (B x + 2 (B 2 + B 3 ) + B 4 ) / 6 (3-40) 

T n+ i = T n ♦ M (C x + 2 (C 2 ♦ C 3 ) ♦ C 4 ) / 6 (3-41) 

B 1 = h-G(T) 

C x = h • F(T ,T) 

B 2 = h G(T+i C x ) 

C 2 = h F(T+i B x , T + I C x ) (3-42) 

B 3 = h G(T*I C 2 ) 

C 3 = h F(T*I B 2 , T + | C 2 ) 

B 4 = h G (T + C 3 ) 

C 4 = h F(T+B 3 , T+C 3 ) 
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This comparison is shown in Table 3-5 for different 
step sizes. 


Numerical Solutions of the Heat Flux and 
Momentum Equations j 

The solution of the heat flux and momentum 
equations are obtained after the temperature profile 
has been established, i.e., the energy equation has 
been solved. The momentum equation was uncoupled from 
the energy equation, as previously explained, because 
the energy dissipated by PV work was considered small 
when compared to the energies by convection and 
chemical reactions. The equations for the heat flux and 
pressure are first order, integral equations with 
variable coefficient. 

The heat flux equation is: 


q NET 




X i 


dT + 



dT 


V H i^i 
i = 1 d T / d Z 


dT 


(3-43) 


The first term in the equation is the heat 
absorbed by convection in the combined decomposition- 
char zone. The second term is that absorbed by the 
polymer degradation in the decomposition zone. Note 
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TABLE 3-5. Comparison of the Fourth Order Runge-Kutta 


Solution With the Analytical Solution of 
2 

— j = A ^ for 100 Integration Steps. 


1 

Numerical 


Analytical 

Distance 

T of 


T of 

Z 

500.0 


500.0 

0 

920.7 


920.7 

0.1 

1385.6 


1385.6 

0.2 

1845.7 


1845.7 

0.3 

1899.4 


1899.4 

0.4 

3094.8 


3094.9 

0.5 

3788.4 


3788.5 

0.6 

4554.9 


4556.0 

0.7 

5402.0 


5402.1 

0.8 

6338.2 


6338.4 

0.9 

7372.9 


7373.1 

1.0 



T = 500°F 

initial 




dT 

dZ = 1 000 

Z=0 
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that the integration limits are from T 0 to T p . Where Tp 
is the temperature at which all of the virgin material 
has been degraded to gases and char. The last term 
accounts for the heat absorbed by chemical^ react ions . 

The pressure equation is: 


, 2 1 

P = |P L + 2R [ J(u/y) W (T/M W ) dZ 

z 


L . 

+ S 6 (T/M w ) W 2 dZ] 

Z 

Equation (3-44) was integrated using Simpson's Rule. 
The general formula for the Simpson's Rule analysis 


(3-44) 


is 


(23) : 


id! = | [f 0 + 4 (f i + £3 + . . . * f 2 n-l ) 

♦ 2 (fj ♦ f 4 . . . ♦ f 2n _ 2 ♦ f 2n ] - 5j£ f* (3-45) 

where 2 ^- f4 i s the truncation error. 

In terms of the pressure distribution equation, 
the Simpson's Rule functions are: 

jf Pi dZ = J*(u/y) W (T/M w ) dZ (3-46) 


Jfp 2 dZ = Jb (T/M w ) dZ 


(3-47) 
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A study on the value of the step size h that 
would minimize the truncation error and maximize the 
accuracy of the approximate solution was done by April 
(2) . A comparison of his results for various step 
sizes is given in Table 3-6. 

Summary of the Theoretical Development of the 
Equations of Change for Flow in tne Combined 
Decomposition Char Zone 

The importance of considering both the 
decomposition zone and char layer as important heat 
absorbing regions was presented. The Equations of 
change (continuity momentum and energy) were developed 
to describe the flow of gases and to predict the energy 
transfer in the combined regions. These equations 
were developed for three flow conditions. Namely 
frozen, equilibrium, and non-equilibrium flow. 

The method of determining the effects of the latter 
two flow conditions was presented. In addition, the 
equations to describe the degradation kinetics of the 
virgin material were discussed. Finally, the numerical 
procedure was discussed, along with the comparison to 
verify the accuracy of the computer implemented solution. 

The following two chapters. Chapter IV and V 
will discuss the methods used to calculate chemical 
equilibrium and non-equilibrium chemistry, respectively. 



TABLE 3-6 Comparison of various Simpson’s rule increment sizes for 
the frozen flow, variable properties model. (2). 
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-Although the development of the chemical equilibrium 
analysis (Chapter IV) will appear to be more difficult 
than for non-equilibrium analysis (Chapter V), the 
computer implemented solution of the latter has given 
severe numerical difficulties which will be thoroughly 
discussed . 
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CHAPTER IV 


CHEMICAL EQUILIBRIUM ANALYSIS 


Introduction 

A general method for computing complex chemical 
equilibrium in a multicomponent, polyphase system is 
presented in detail (1, 2). This method, called the 
Chemical Equilibrium Analysis, is a modified version of 
the RAND method, which is also known as free energy 
minimization. It was originally presented by White, 
Johnson and Dantzig (3) for computing complex chemical 
equilibrium for an ideal gas system. It was subsequently 
extended to multiphase systems by Hubert and Stephanou 
(4) and others (5-7) . In addition to this discussion 
of the method, brief comparisons with some other well 
known methods (8, 9, 10) for computing chemical 
equilibrium are presented. These comparisons will 
justify the choice of free energy minimization as a 
more flexible and easier method to implement on a 
machine for highly complex chemical systems. 

The literature on the computation of complex 
chemical equilibrium is so extensive that it would 
detract from the main purpose of this chapter to review 
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all of it especially since a very extensive and well 
documented review on the historical development of 
chemical equilibrium was published recently by Zeleznik 
and Gordon (11, 12). A more modest undertaking was 
presented by del Valle et al . (13) . At the end of 
the chapter, comments will be made on some of the other 
publications (16-66) which are available for solving 
the nonlinear chemical equilibrium equation. It will 
be realized that all these techniques are but minor 
variations of the most general methods (3, 8, 9, 10). 

It should be noted that the solution of the chemical 
equilibrium problems can be categorized either as 
the solution of a system of non-linear equation by 
the functional iteration mpthod, i.e., equilibrium 
constant formulation (8, 9) or, as the direct 
minimization of the Gibbs free energy by descent 
method (3, 10). In this chapter the latter method is 
developed . 

Condition for Chemical Equilibrium 

The condition for chemical equilibrium in any 
reacting system at constant temperature and pressure 
is that the change in free energy of the reactions 
must be zero (14). Stated another way, the sum of 
the free energies of the individual components, i.e. 
the total free energy of the system, shall be a 
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minimum. A negative free energy may be considered 
as a driving force which causes a reaction to approach 
equilibrium. It can also be considered to be a 
measure of the departure of the reacting system from 
an equilibrium state. Thus, it is through this 
departure and hence, through its minimization that 
chemical equilibrium is calculated. 

Stochiometric Balance 

The law of the conservation of mass requires 
that for a chemically reacting system with no nuclear 
reactions that the chemical elements be conserved 
regardless of the number of phases being present at 
equilibrium. The stochiometric relation can be 
expressed as: 

n 1 

T a'ij x i + Z a ij xi = bj j = 1 m (4-1) 

i= 1 i=n+l 

where a^j is the formula number giving the amount of gram- 
atoms of the jth chemical element in the i t ^ species. 

The moles of each specie present in the system is denoted 
by xi. The gram-atoms of each element j is denoted by 
bj . The number of gaseous chemical species is n, and 
the number of condensed species is denoted by n+1 to j . 
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Conservation of Charge 

When ionization is present in a system, the 
conservation of charge can be expressed in a form 
similar to Equation (4-1). Thus: 


1 


l 

i = l 



(4-2) 


where a^ is the charge on the i*^ specie, is 

’m+l 

the number of moles of each specie, and is zero, if 

the system is electrically neutral. Thus if we 
consider a charge as the m+l element we can use an 
equation of the form of Equation (4-1) . 


The Equation of Free Energy 

It was previously stated that one of the 
conditions for chemical equilibrium was that the 
total free energy of the system be a minimum. To 
minimize the free energy an appropriate expression 
of this function is necessary. At equilibrium, the 
free energy of each specie in a mixture at any 
temperature T can be obtained by integrating the 
definition of fugacity (14) in terms of the free 
energy (dFp = RTdlnfi) . This integration gives: 


( F T) i - (F^)i = RT In Ti/^ 0 i = i... n (4-3) 


but 
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Equation (4-4) is the definition of activity. Thus, 
Equation (4-3) becomes: 

» 

(F T )i - (Ft 0 ) i = RT In aj 

or i = 1 , . . .n (4-5) 

(Ft) i = (F T 0 )i + RT In ai 

lip to this point there have been no restrictions 
posed on a chemically reacting system. At this 
juncture the two restrictions are introduced. One is 
the restriction that the chemically reacting gases 
are ideal. This in itself is not too severe if the 
analysis is constrained to high temperatures and low 
pressures, which is precisely the area of investigation 
in this research.* The second restriction is to 
assume that there is no mixing among the condensed or 
solid species, i.e., they remain in their pure state. 
Hence the activity of the solid can be taken to be one. 

*In references 6, 45 and 46, non-idealities 
are considered. 
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Free Energy Minimization 

Ideal Gas . For an ideal gas ct^ the activity of 
the specie i can be taken to be equal to its partial 
pressure, pi. Thus Equation (4-5) becomes: 

(Ft) i = (FT°)i + RT In pi i = l,...n (4-6) 

Solid Phase . In the event of solid carbon or 
any other condensed specie present in the system, 
Equation (4-5) can be simplified by recalling the fact 
that the activity of a solid or pure condensed specie 
is one. Therefore: 

(Ft) i = (F T °)i i = n+1 , . . . 1 (4-7) 

Equations for Free Energy of a Complex Gas - 
Solid Mixture" . Consider admixture containing 1 chemical 
species (n gaseous and n + 1 . to 1 condensed) from 
chemical elements. Let 

x = l \ i (4-8) 

i=l ’gas 

where is the moles of specie i of the gas and x is 
the total moles of the reacting mixture. The definition 
of partial pressure is: 

Pi = P . i = 1. . .n (4-9) 

x 
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where P is the total pressure of the reacting system. 
Substituting Equation (4-9) in (4-6) and multiplying by 
xi/RT results in: 


xi (Fp) , xi (F t 0 ) j _ 

— Rf = — RT + Xi t ln P + ln x i / x 1 


i = 1 . . .n 


(4-10) 


It is convenient to define 


and 


f. = x i ( p T) i 

L i 

RT 


i = 1 


.1 


( FT °) 1 + ln P i = 1 , . . .n 

RT 


(4-11) 


(4-12) 


Combining Equations (4-11) and (4-12) with (4-10) gives 
the free energy function for the gases and is shown 
below: 

f i , = xi (ci + ln i = l,...n (4-13) 

gas x 

It can be observed that ci will be constant at constant 
T and P, and xi will be the only variable to be 
determined from thermodynamic equilibrium. A similar 
expression for the solid or, condensed^ phase is obtained 
by multiplying both sides of Equation (4-7) by 
xi/RT. Making use of the definition of Equation (4-11) 
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results in the free energy function for the solid, which 
is 


f 


i , solid 


Xi (F T °/RT) i 


i = n+1 , . . .1 (4-14) 


The total free energy of the mixture is obtained by 
summing Equations (4-13) and (4-14), and this gives: 

n 1 


F(X) = l 

i = l 



1 fi>solid 
i=n+l 


(4-15) 


or 


n 1 

F (X) = l xi [ci + ln xi/x] + l xi(Fx°/RT) i 
i=l i=n+l 


(4-16) 


where X = (xi , x2,...,x ) 

A set of equations is now developed to enable the 
computation of the set of xi's that satisfy Equation 
(4-1) (material balance) and makes the free energy 
function, of Equation (4-16) a minimum. 

Rationale for the Iterative Procedure . At a 
given temperature and pressure it is necessary to 
determine the amount of each chemical specie, xi , 
present that minimizes the free energy. However, 
problems arise when standard minimization techniques 
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are applied directly to F(X) and the associated 
material balance constraint. This is because the 
equations resulting from setting the first partial 
derivation with respect to xi equal to zero cannot be 
solved explicitly, for the independent variables, xi. 
Therefore, it is convenient to make a quadratic 
approximation to the free energy function which 
permits the independent variables to be expressed 
explicitly. Then an iteration scheme is developed 
such that this quadratic approximation approaches 
the actual value of the free energy function at the 
point of minimum free energy. To develop these 
equations for the quadratic approximation and the 
iterative scheme, select any positive set of mole 
values of Y = (yi, y2»«**yl) which satisfy the 

material - baTance— equations-^— Equation— (-4— -1-)-,- -as— ' the 

initial estimate of the composition. There, the 
value of the total free energy function of the 
mixture is: 

F (Y) = l yi (ci + In + I Yi (F T ° ) i/RT (4-17) 
i=l > i=n+l 

However, it is not necessarily true that 
the assumed mole numbers are the ones for the 
components in thermodynamic equilibrium. For this 
to be true the free energy must be a minimum and the 
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material balance equations satisfied. Therefore the 
objective of the iteration scheme to be developed 
is to locate the point of minimum free energy by 
using the quadratic approximation to the free energy 
function and have the material balance constraints 
satisfied. The approach will be to form the quadratic 
approximation of the free energy function at X 
expanded about Y. Then the adjoint equation 
employing Lagrange multipliers will be formed with 
this approximation and the m material balance con- 
strait equations. At this point, the usual procedure 
would be to set partial derivatives with respect 
to the independent variables, (xj_, X 2 ,*».x^) and the 
Lagrange multipliers ( tt ^ , tf2»‘** 1T m) equal to zero, 
and solve simultaneously the resulting set of m + 1 
equations for the minimum of the constrained quadratic 
approximation. Then the procedure would be repeated 
employing some convergence criteria to approach the 
point of minimum free energy. 

It is more convenient, however, to use a 
variation of this procedure which requires that only 
m + 1 + s equations be solved simultaneously, where 
s is the number of solid or condensed species, rather 
than m + 1 + n + s equations. Since a quadratic 
approximation is used, the linear equations resulting 



143 


from taking the partial derivatives with respect 

to the xi ' s can be solved directly for the xi'| as in 

terms of the Lagrange Multipliers, tt . *s and x, 

3 

the total moles of gases. This permits eliminating 

x i gas *he material balance equations and results 

in the m + 1 + s equations to be solved simultaneously 

for the it *s x, and x^, solid. There are m equations 
3 

with the ffj ’s as the independent variables which 
result from the elimination of xi in the material 
balance equations on the individual elements, Equation 
(4-1) . The other equations are Equation (4-8) and the 
ones from setting the partial derivative with respect 
to x i, solid of the adjoint equation equal to zero. 

These equations will be developed and placed into a 
matrix form. Then a convergence procedure will be 
discussed to insure that the point of minimum free 
energy is obtained. 


Quadratic Approximation to the Free Energy 
Function. The problem is to expand F(X) in the 
neighborhood of Y in terms of the difference 
A i = ( x i " y^) • T° do this a Taylor series 
expansion for the free energy function at X, F(X), 
expanded about Y, the estimated values, results in 
the quadratic approximation to F(X), which follows: 
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+ l . Hf Y) * A i + 7 l , l , ( Ai + Ak 


Q(X)"F(Y) - j ayj 


i=l k=l 


3yl * ^7k 


F (Y) + 


(4-18) 


where 


( A i 3 y i + Ak 3yk ^ 


A- 2 3 2 -, + 2 A- A, — 

1 sy i 2 1 k 3yi 3 yk 


2 3 


k 3yk 1 2 


(4-19) 


Neglecting any terms higher than second order. Equation 
(4-18) becomes: 


1 1 2, 


1 ”3F (Y) - . if ? ' 3 ^ F( Y) A ; A .- 

<X X) = F (Y) + J = i ly— A a ♦ - l = i I = i ayTay^ A i A k 


(4-20) 


Taking the appropriate partial derivative gives: 


|£Ci= Ci M. yi/y 

3y x 


1 < i < n 


(4-21) 


3F(Y) _ (F T °)i 
3y i ' RT 


n ♦ 1 < it 1 


(4-22) 


3 2 F(Y) _ 6ik . 1 

3yi Yk yi y 


1 < i < n 


(4-23) 
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where 6^ is the Kronecker delta (6^. = 0 for i = j 
and 6 ilc = 1 for i^j) and not that the second partials 
in the case of condensed species are zero. 

Substituting Equation (4-21) -, (4-22) and 
(4-23) in (4-20) results in the quadratic approxima- 
tion to the free energy function at X expanded about 
Y, and is: 


Q (X) = F(Y) + 


n 


l 

i = l 



+ In + 

y 



* I! I (^-i) Ai Ak 
z i=i k=i yi y 


(4-24) 


where Q(X) is the quadratic approximation to the free 
energy function. It has been shown that both F(X) 
and Q(X) are positive definite and hence, both 
functions are convex (3, 4, 11, 53). Thus, each 
has a unique minimum. 

Lagrange Multiplier Formulation and Minimization . 
The problem is now one of obtaining the minimum of the 
function Q(X) subject to the material balance con- 
straints of Equations (4-1). The adjoint equation, or 
constrained equation, is formed using the method of 

Lagrange's undetermined multipliers. The result is: 

m n 1 

G(X)=Q(X) + l tt, (bj - I aij xi - l aij Xi ) (4-25) 
j=l i = 1 i=n+l 

where irj is the Lagrange multiplier. 
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At this point the usual procedure would be to 
set partial derivatives with respect to the independent 
variable (xi , x 2 , x 3*** x l) the Lagrange multipliers 
c rr 1 , ^2 , . . . u m ) equal to zero. That is, to minimize 
G(X) it is required that the following conditions 
be satisfied: 


3G 

9x j 


0 


and 



(4-26) 


Then the procedure would be to solve simultaneously 
the resulting set of m+1 equations where 1 = n + s, 
s being the total number of solid or condensed species. 
This would be. repeated employing some convergence 
criteria to approach the point minimum free energy. 
However, as has been said before, it is more convenient 
to use a different procedure which requires that only 
m + 1 + s equations be solved simultaneously rather 
than m + n + s equations. Since a quadratic 
approximation is used, the linear equations resulting 
from taking the partial derivatives with respect to 
the x^'s can be solved directly. Performing the 
partial differentiation of Equation (4-25) ; the 
constrained quadratic approximation, results in the 
following expression: 
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mm 

i£OQ = 3 Q( X ) - l l aij n j = 0 1 < i < 1 ('4-27) 

3x i 3x i j=l j=l 

3Q (X) 

It is necessary to develop the expression for 

Thus, taking the derivative of Equation (4-24) results 

in: 

= ( ci + In yi/y) + ( ~ - ) 

Xi Y i .7 

1 < i < n (4-28) 

i = n + 1 1 (4-29) 

Substituting Equations (4-28) and (4-29) in (4-27) 
results in the following expressions: 

^ r* 

_ = [ Ci + In yi/y] + [xf/yi - x/y] 

m 

= 0 i = 1 n (4-30) 

m 

- J tt j aij = 0 i = n + 1...1 (4-31) 

j = l 


S j "j a u 


and 


3G = (F°x) i 
3x i RT 


and 

mil « Ft 0 

3x t R T~ . 

x l 


Since yi > 0 we may solve for xi in (4-30) and obtain: 
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x i ,gas 


m 

fi(y) + (yi/y) x + l (it* a i: j) Yi 

j=l 3 


i = 1, . . .n (4-32) 

This gives n explicit equations for the moles of the 
gases in terms of the m + 1 unknowns, and x. Equation 
(4-23) permits the elimination of x^gas from the 
material balance equations. Substitution of Equation 
(4-32) into (4-1) results in m equations of the form 
given below: 


n 


I 

i = l 


ajk yi 

y 


x 


m n 

+ l I [Uij) ( a ik) yil 71 j + 
j-1 i-l 


1 n 

l a ik xi = bk + l aik fi (Y) k = 1 

i=n+ 1 i=l 


,m 


(4-33) 


For the purpose of simplifying the above, it is 
convenient to define: 


r ' 

l a ik k = l,...m (4-34) 

and 

x 

y 


(4-35) 
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In addition 


n 

r jk = Tkj = ( a i j) (aik) Y i (j> k = 1 > 2 > . . .m) 

(4-36) 


Equation (4-33) on re-arranging becomes: 


I rjk tt j + b k u + l a ik x i = b k + 
j = l J i=n+ 1 


l a ik f, (Y) (k * l,...m) (4-37) 

i = l 

This represents m equations in m + 1 + s unknowns 
irjs, u and so lids* An additional independent 
equation is obtained by noting that = x and by 

summing Equation (4-32) from i = 1 to n the following 
equation is obtained. 


m ( n 

l b- tt = l f L (Y) (4-38) 

j=l J J i=l 


The additional s equations needed to compute the 
amount of condensed (or solid) phases present are 
obtained from Equation (4-31). This results in: 


m 


(F T °/RT)i = l tt, 

j = l J J 


i = n + 1. . .1 (4-39) 
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These equations, which are to be solved simultaneously, 
are summarized and given in general form in Table 4-1. 
When the equations are solved, values for the moles 
of condensed species are obtained directly. Also, 
the values of the Lagrange multipliers ttj and that 
of u are obtained. These two latter values are 
substituted in Equation (4-32) to calculate the 
moles of the gas species. A more detailed exposition 
of the procedure to follow is given below. 

Computation Procedure . The iterative pro- 
cedure is initiated by assuming any positive solution 
Y - (yi, y2>***yi) which satisfies the material 
balance equations, Equation (4-1). Then, the value 
of fi(y) are determined by Equations (4-13) and (4-14) 
as are the “values ~of r by-Equa-tion -( 4- 36) .. Next 
the system of equations in Table 4-1 are solved 
simultaneously. Note that the solution to these 
equations give m values of the Lagrange multipliers. 
tt j , the value of u = (x/y) and values of x^'s solid. 

In addition, it should be noted that the values of 
x i, solids are t * ie onl y ones directly obtained from 
the simultaneous solution of these equations. The 
values of xj^gas are not obtained directly from this 
solution, but must be computed using Equation (4-32) 
which is in terms of the known values ttj and u. The 




TABLE 4-1. General Equations for the Solution of the Equilibrium Composition of 
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n+ljl^i + • • • + a n+q,m 7T m : (F°/RT)k+1 



TABLE 4-1 (continued) 
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new and improved values of can be used, if they are 
all positive, as a starting point for the next 
iteration. However, negative xi's can occur and steps 
should be taken to avoid this possibility and thus, 
insure convergence. This is discussed below in the 
context of common numerical problems encountered 
when using iterative schemes. 

Numerical Problems in Convergence Procedure 

Convergence in an iterative calculation, 
involves, usually, three numerical problems: (1) How 

to insure numerical convergence. (2) How to determine 
the step size for a given iteration. (3) How to 
determine when to stop in the computation cycle. 

How, to Insure Numerical Convergence . Normally 

in the iterative procedure, the amount of each specie 
x^ which is calculated at the minimum of the con- 
strained quadratic approximation, is used as the new 
estimate for the following iteration. To insure that 
oscilations and over-corrections will not occur, the 
following convergence scheme similar to that of White 
et al . (3) was used. 

To prevent oscillations, the gradient of the 
free energy function can be examined, and the values 
of the x^’s determined that insure the free energy 
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will always decrease. To do this, a parametric 
representation of the free function is formulated. This 
is done by taking the values of the moles of each specie 
calculated at a given iteration to estimate the values 
of the moles of the chemical species yi( ne w) f° r t ^ ie 
next iteration. Thus, 


Y i (new) ~ ^i (old) + 


1 < i < 1 (4-38) 


and 


A i “ x i ^i(old) 


1 < i < 1 (4-39) 


where A is the parameter of the line through xi, and 
^i,old» an ^ it can vary from zero to one. The value 
for 1 that was used to insure convergence is discussed 
in _a later part of this section. 

The free energy function can be expressed in 
terms of the convergence parameter A by substituting 
yi (new) of Equation (4-38) in Equation (4-16). The 
result is: 


n 


FU) = l [yi + AA-J [Ci + In yi + **Aj/y +Ad ] + 
i = l 


L, lyi * XAi! • (m)i 


(4-40) 


where in the above expression yi is equivalent to yi(old) 
of Equation (4-39) and A = x - y. 
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To determine the direction (increase or decrease) 
of the free energy function, the derivative with respect 
to X is examined after every iteration. The derivative 
is easily computed and is: 


dF(X) 

dx 


n 


1 ^ [ Ci + i n 

i = l 


yi + *Ai /y + X aj 


+ 


1 


l 

i=n+ 





(4-41) 

If the directional directive dF/dX is negative 
a descent path is followed. That is the value of the 
free energy function in the r+1 iteration will be less 
than in the r*h iteration. When this procedure is 
continued the minimum is eventually reached. However, 
when a non-negative value of the directional derivative 
is obtained, the value of X (step size) should be 
reduced until a negative value of this derivative is 
obtained. A more thorough discussion on the proper 
value of X is given below. 

How to Determine the Step Size for a Given 
I teration . The original paper of White, Johnson and 
Dantzig (3) on which the technique developed in this 
chapter is based, does not shed much light on the size 
of the iteration step, or convergence parameter X that 
is required for converences . However, recently Zeleznik 
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and Gordon (11, 12) showed that for a function of A 
F(A) s [F (y + X A) ] * the optimum value of the step 
size can be found. This is possible since the condition 
for equilibrium is that the free energy be a minimum. 
This is equivalent to determining X by solving the 
equation : 



y n ew = ^oid * x A i 

where &T denotes the transposer of the matrix. 

This optimum value of X allows for a better 
control of the path of descent with the corresponding 
decrease in the number of iteration necessary to reach 
the minimum. 

Practical consideration, however, dictate that 
a minimum amount of time be spent in this phase of the 
iterative procedure. Consequently, it is best to 
sacrifice a few extra iterations and obtain in return 
a faster estimate of the value of X. This was done 
by Oliver et al . (15) . Their technique consisted in 
choosing the specie with the largest negative value of 
hi and computing a value of x' such that yi(new) be 
zero; therefore 


* The bar ( ) denotes a vector. 
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• _ . yi(old) (4-43) 

- 

Then X was set to some fraction of X*. They proposed 
0.99 X*. We have found from experience that 0.9X 
works better. 

However, this procedure alone, does not insure 
convergence. It is possible to fine a value of X 
which meets the above requirement but yet makes 
dF/dX positive. This has the effect of increasing the 
free energy rather than decreasing it and can cause the 
solution to diverge. Therefore, one must check if X' 
calculated by Equation (4-43) makes the gradient 
dF/dX negative. If it does not, the value of X' 
must be further reduced until this condition is 
satisfied. 

It is possible to encounter situations where 
the value of X becomes practically zero. This will 
most likely occur when trace species are present, that 
is, when the value of yi , the moles of a specie is very 
small. A criteria has been developed to automatically 
eliminate trace species from the computation in the 
chemical equilibrium analysis program. These procedures 
are further discussed in the section on trace species. 

It should be stated that one should never allow the 
correction XAi to be less than 6 or 7 orders of 
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magnitude smaller than y^ since this would make 'the 
convergence scheme very time consuming. More will also 
be said about this in the section dealing with trace 
species . 

It has been observed that as the minimum is 
approached large changes in A occur, which cause the 
value of dF/dA to become smaller. This allows for a 
control path of descent, and convergence to the minimum 
is obtained. 

How to determine where to stop . It was 
established in the above paragraph that as long as dF/dA 
remains negative, the iterative cycle will continue to 
calculate new and improved values of the y's. In 
principle the iteration should be stopped when: 

" n - - - - 

I I Ai | = 0 (4-44) 

i=l 

However, this criteria is impractical to achieve 
in a digital machine because of round-off error. 

Rather, the criteria should be determined by these 
four basic questions: 1) How accurate should the 

answer be? 2) How many significant figures are available 
in the digital machine? 3) Can a little accuracy be 
sacrificed for time? 4) How accurate is the input data? 
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most cases the data is the constraining 
factor. There is no point in getting successive 
checks in the six or seven significant figure when the 
thermodynamic data is only good to three at most four 
significant digits. We have found from experience 
that when 

1 

0.001 < l | A ± | < 0.01 (4-45) 

i 

rapid convergence and accurate values result. 

There are other criterias that can be used to 
test for convergence. One would be to take successive 
values of the free energies and compare them. When the 
difference is less than a pre-requisite number of 
significant figures stop the iteration. Similarly, 
the value of the gradient dF/dX can be -examined-.- - — - — 
When its value is close to zero we are assured of 
having reached the minimum, i.e., equilibrium. 

Although these last two criterias seem to be logical 
choices, we have found that when using Equation (4-45) 
better convergence control is achieved. 

Gaseous Trace Species 

When very little is known about a particular 
chemical system, one is forced to consider as many 
chemical species as thermodynamic data is available 
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for them. In the case, for example, where carbon and 
hydrogen react there are many possible products of 
reaction. At high temperature it would be expected 
that lower molecular weight hydrocarbons would exist 
while at lower temperatures, higher molecular weight 
species would be expected. This is probably the 
most one can safely infer from basic chemistry know- 
ledge. But little else can be said. Hence, when it 
is desired to know the equilibrium composition of a 
hydrocarbon mixture at different temperature and 
pressure, one, of necessity, needs to consider all 
probable hydrocarbon species. Under these conditions, 
it is very possible that many of the chemical species 
assumed to be present will only be in trace quantities. 
However, it will be illustrated how species which are 
trace at a particular temperature aifd~ pressure can 
become important at other temperatures. This will show 
that trace species cannot arbitrarily be eliminated 
from the computation, especially, for those systems - 
where varying temperatures and pressure can be en- 
countered. Thus, it becomes important to devise an 
automatic procedure to eliminate trace species momentarily 
from the computations so as to speed up the convergence, 
to the minimum free energy, and later reintroduce 
them at the end of the iteration. 
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To show that the elimination of trace species 
has no effect on the equilibrium composition of major 
components. Table 4-2 illustrates this with the example 
of the equilibrium composition of a 40 percent nylon 
and 60 percent phenolic resin composite. The tempera- 
ture is 800°C and the pressure one atmosphere. Run 1 
shows the equilibrium composition for twenty one assumed 
species. Run 2 was performed with C 2 H, C 3 II, C 4 H', and 
C 3 taken as zero composition. Since the equilibrium 
composition of Run 1 predicts that these four species 
have negligible mass fractions, their elimination 
produced no effect on the equilibrium composition 
of the seventeen remaining species. 

In Table 4-3, however, the results are 
strikingly different between Run 1 and Run 2 for a 
temperature of 3100°C. In this case the elimination 
of the species C 2 H, C 3 H, C 4 H and C 3 does have a 
profound effect on the composition. Although these 
four species were trace species at 800°C this is now 
not true at 3100°C, and thus illustrates the importance 
of considering as many species as thermodynamic data 
is available. This is especially true when the 
chemical system under consideration is studied under 
a diverse range of temperature conditions. 
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TABLE 4-2. Equilibrium Composition of the Pyrolysis 
Gases of a 40 Percent by Weight of Nylon 
and 60 Percent Phenolic Resin at P = 1 atm 
and T » 800°C. 


Mass Fraction 

Species Run 1 Run 2 


H 

0.5483 

X 

10' 9 

0.5483 

X 

10" 9 


0.6788 

X 

10' 1 

0.6788 

X 

10*1 

CM 3 

0.1148 

X 

10’ 7 

0.1148 

X 

10* 7 

ch 4 

0.1893 

X 

10 _1 

0.1893 

X 

10" 1 

c 2 h 2 

0.7475 

X 

10' 8 

0.7475 

X 

10*8 

C 2 »4 

0.5910 

X 

10' 6 

0.5910 

X 

10" 6 

C 2 H 6 

0.4606 

X 

10‘ 6 

0.4606 

X 

10" 6 

C 6 H 6 

0.8471 

X 

10-13 

0.8471 

X 

10-13 

n 2 

0.3875 

X 

10 " i 

0.3875 

X 

lO'l 

Nil 

0.3275 

X 

10" 4 

0.3275 

X 

10' 4 

HCN 

0.4101 

X 

ib* 5 ~ 

" ' 0.4101 

X 

10- 5 

h 2 o 

0.1409 

X 

10 " 1 

0.1409 

X 

io" 1 

OH 

0.8392 

X 

10-12 

0.8392 

X 

10*12 

C02 

0.8776 

X 

10-2 

0.8776 

X 

10 " 2 

CO 

0.2220 



0.2220 



CN 

0.1006 

X 

10-14 

0.1006 

X 

10" 4 

c 2 h 

0.1347 

X 

10-16 

0 . 



c 3 h 

0.6904 

X 

10-18 

0 . 



c 4 h 

0.5722 

X 

10-22 

0 . 



C3 

0.7290 

X 

10-28 

0 . 



C 

0.6295 



0.6295 
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TABLE 4-3. Equilibrium Composition of the Pyrolysis 
Gases of a 40 Percent by Weight of Nylon 
and 60 Percent Phenolic Resin at P = 1 atm 
and T = 3100°C. 


Species 

Run 

Mass 

1 

Fraction 

Run 

2 

H 

0.2510 



0.2859 



«2 

0.3318 



0.4304 



ch 3 

0.9268 

X 

10- 4 

0 .1369 

X 

10' 3 

CII4 

0.5311 

X 

10‘ 5 

0.8937 

X 

10' 5 

C 2 H 2 

0.6010 

X 

10' 1 

0.7796 

X 

10 _1 

c 2 m 4 

0.1526 

X 

10‘ 5 

0.2567 

X 

10' 5 

C2H6 

0.7115 

X 

10- 11 

0.1553 

X 

10-10 

C6H6 

0.4456 

X 

10' 21 

0.9725 

X 

10‘ 21 

N2 

0.3377 

X 

10’ 2 

0.8633 

X 

10* 2 

NII 3 

0.9244 

X 

10-7 

0.1386 

X 

10‘ 6 

HCN 

0.2175 

X 

10' 1 

0.2514 

X 

10-1 

h 2 o 

0.2560 

X 

10' 6 

0.3611 

X 

10~ 6 

OH 

0.8525 

X 

10-7 

0.1056 

X 

10~ 6 

C0 2 

0.1467 

X 

10-7 

0.1735 

X 

10' 7 

CO 

0.1511 



0.1643 



CN 

0.7421 

X 

10" 2 

0.7533 

X 

10' 2 

c 2 h 

0.7364 

X 

10' 1 

0 . 



c 3 h 

0.5525 

X 

10 -1 

0 . 



C4H 

0.3552 

X 

10" 1 

0 . 



C3 

0.3938 

X 

10’ 2 

0 . 




C 


0.2562 


0.7603 
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It is also important to note that the speed of 
convergence can be drastically affected by whether or 
not trace species are eliminated during the iterative 
procedure to minimize the free energy. For the 
example in Table 4-2 if all of the species with mass 
fraction less thatn 10"^ are eliminated the number of 
iterations needed to converge to the minimum were 
seven, compared to twenty-one when the twenty-one 
species were considered. 

It is possible to eliminate trace species during 
the convergence scheme, and when convergence has been 
achieved, recompute the value of trace species. This 
can be calculated by using the following equation: 

m 

= x exp (-Ci + l a-s tt ^ ) i » 1 . . .n (4-46) 

j = 1 

Equation (4-46) is obtained by differentiating 
Equation (4-24) with respect to x^ and replacing Q(X), 
which is the quadratic approximation of the free energy 
function, with the exact function F(X) . Note that the 
total number of moles x and the Lagrange multipliers 
Vj are constants since trace species have no effect 
on these values. Hence, the only necessary thing to 
compute xi is to compute ci from Equation (4-12). 
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Solid or Condensed Phase Species 

It has been our experience that apart from 
control of the minimum size of the value of the con- 
vergence parameter X, there has been no numerical 
difficulties encountered in determining equilibrium 
compositions in gaseous systems when trace species 
are eliminated during the iterative procedure. As has 
been previously stated, these trace species are 
reintroduced at the end of the iteration and their 
numerical values are recomputed by the use of Equation 
4-46. Straight forward application of the iterative 
equations in Table 4-1 results in a rapid convergence 
to the minimum. In principle the same applies to 
systems with pure condensed species. However, 
difficulties are known to occur in these computations, 
and care should be exercised to circumvent them as 
will be discussed. 

When formulating the chemical equilibrium 
problem for a particular chemical system, one assumes 
a number of chemical species that might be present at 
equilibrium. Some of these species will turn out to 
be present in trace quantities. Equation (4-46) is 
used to compute the numerical values of the moles of 
each of these trace species. Unfortunately, Equation 
(4-46) is not applicable to other than gaseous species. 
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Consequently, this procedure cannot be used for solid 
or condensed species. There is, however, a procedure 
and certain rules derived from thermodynamics that 
need to be followed rigorously for the elimination or 
addition of condensed species. These procedures are 
explained below. 

When to Eliminate a Solid Species . There can 
be times during the convergence scheme when the moles 
of a solid or a condensed phase species becomes very 
small or negative. When this occurs the value of the 
convergence parameter X, needs to be reduced. If 
this condition persists the value of X can become 
zero. When this occurs the iterative procedure might 
not converge. It then becomes necessary to eliminate 
this condensed specie to assure convergence. 

When to Add a Solid Specie . The complimentary 
problem of adding a condensed specie must be approached 
from a different viewpoint (4) . Consider the adjoint 
Equation (4-24) . This is the expression for the total 
free energy function to be minimized. Since Q(X) = F(X) 
at the minimum Equation (4-24) becomes: 

m n 1 

G (X) = F(X) + l TTj (b. - l aijxi - l 

j=l J i=l i=n+l 


aijxi) (4-47) 
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when the material balance constraint is satisfied: 


G (X) = F(X) (4-48) 

Suppose that for a given set of mole numbers, 

t » t f 

X = (x-^, x 2 . . • x^) , F (X) is a minimum. Then: 

i = 1 . . . n (4-49) 

i = n + 1 . . .1 (4-50) 

Now consider adding a condensed specie to a 
system where this specie is already present in the gas 
phase. Denote this as the q' condensed specie and x' 
as the number of moles of the specie. The change in 

free— energy •- due- to the- addition, of this ..spe.c.ie__to 

the system can be obtained using a Taylor's series 
-expansion, and is: 


3G_ 

3X| = 0 

X=x ' 




n 


l 

i = l 



+ 


1 

I 


i 


=n+ 


1 


3 G 
3x 


a : 

i 


where 


a; 


= Xi 


x: 


(4-51) 


l 


(4-52) 
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Combining Equations (4-49), (4-50), (4-51), 

gives : 

a n _ 3G , , 

AG " JZ* A q (4-53) 

For the free energy of the mixture to decrease, 
3G/3x' must be negative. To minimize the free energy, 
the partial derivatives of the free energy function with 
respect to mole numbers are set equal to zero. For 
the condensed specie this gives: 


H~ = [F ? /RT] i ' 


= o 


(4-54) 


Therefore, if the added condensed specie is to 
decrease, the free energy must have a value such that: 


m 


l ^ 7t j a A j > (Tr/RT)i 


(4-55) 


If the above is true for an added condensed 
specie, then it should be considered in the chemical 
system as the free energy will be less with it 
included . 


Comparison of Equilibrium Composition Calculations 
Witn Published Data and Computations 

In this section comparisons of published results 

with the results of the chemical equilibrium analysis 
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program are presented. Comparisons are made for five 
chemical systems. These are CO, CO 2 , 0 and O 2 system; 
the ammonia equilibrium, rocket fuel reaction product 
system, the me thane -water system and the hydrazine 
oxygen reaction products. 

CO 2 , CO, 0 and O 2 System. A comparison of the 
calculated values of the chemical equilibrium analysis 
program with results published by Kondratev (56) for 
a CO, CO 2 , 0 and O 2 is shown in Figure 4-1. The 
continuous line represents the equilibrium composition 
as calculated by the equilibrium program. The dotted 
lines represent the values given by Kondratev (56) . 

The agreement is excellent. The slight difference is 
due, probably, to the fact that different sources of 
thermodynamic data were used. 

The Ammonia Equilibrium . A comparison is 
given of the equilibrium composition of the N 2 , H 2 and 
NH 3 system as calculated by the equilibrium analysis 
program and the data presented by Dodge and Larson (57) . 

The comparisons are presented for Tables 4-4, 4-5, 

4-6 and 4-7 for pressures of 10, 30, 50 and 100 
atmospheres, respectively. The data presented in these 
tables is the actual experimental values of the equilibrium 
composition of ammonia as measured by several individuals 
and presented by Dodge and Larson (57) . These results 
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1100 1500 1900 2300 2700 3100 3500 3900 

Figure 4-1. A Comparison of a Calculated Values of the Equilibrium Concentration 
Predicted by the Free Energy Minimization With Those of Kondratev (56) 






TABLE 4-4. Equilibrium Composition (Mole Percent) of Ammonia at 10 Atmospheres. 
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TABLE 4-5. Equilibrium Composition 1 (Mole Percent) of Ammonia at 30 Atmospheres. 


172 


E 









aj 

rH 

oo 

to 

r-- 

CT* 

CM 

DO 



rH 


a> 

CH 

CM 

00 



GO 

• 

• 

• 

• 

• 

• 

• 


O 

CD 


o 

oo 

vO 

^1* 

DO 



rH 

rH 

rH 






a, 










<M 




o 


CM 


If) 

CM 


vO 


to 


O 



• 


• 


• 


• 


0 

rH 


rH 


vO 


to 


z 

CM 


rH 







iH ( 


CM 


to 


CM 


4) 

00 




rH 


vO 


rG 

• 


• 


• 


• 


03 

00 


o 


vO 

* 

to 



rH 


rH 






« 









-J 









• 

CO 


LO 


v£5 


CD 


a 

r- 


rH 


oo 



0 

• 

• 


• 


• 


• 

u 

z 



o 


LO 


to 

3 


rH 


rH 





4-> 

a 








rH 









X 3 
u o 

O -H 
4-> t-c 









rt oo 
u < 









o 

J3 4H 

<u 

o 

LO 

CD 

cn 

rH 

to 

oo 

os o 

GO 

oo 

to 

o 

LO 

00 

LO 


a 


• 

• 

• 

• 

• 

• 

• 

+J 

O 


to 

o 


LO 


to 

X c 

a 

rH 

rH 

rH 





u 0 
>- E 
rt 4-» 
0) S-c 
co rt 
0 a 









a 0 

t — N 








a 

u 








c 

o 








<D 10 

' — * 








GO 0 
O -M 

0 








^ as 

Ih 








4-> +j 

3 








•H CO 

+J 

o 

LO 

o 

LO 

o 

LO 

o 

z 

03 

LO 

r*. 

O 

CM 

LO 


o 

T3 

s-< 

to 





■^r 

LO 

03 0 

0 








0 +J 

a 








X -H 

£ 








•H 3 

0 








a 3 

H 








* 



TABLE 4-6. Equilibrium Composition (Mole Percent) of Ammonia at 50 Atmospheres. 
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show that the chemical equilibrium analysis program 

! 

can predict the chemical equilibrium composition for 
systems in equilibrium. It is suspected that part of 
the difference between the actual and calculated 
values are due to non- ideali ties and to the experimental 
error of the measurement. In general the results 
agree within 51 . 

Rocket Fuel Reaction Products . In Table 4-8 
a comparison of the equilibrium composition of unsymetri 
cal dimethyl hydrazine (UDNII) , red fuming nitric acid 
propellant combination as calculated by the equilibrium 
analysis program and by Brandmeir and Harnett (32) 
is presented. 

These authors calculated the equilibrium 
composition using free energy minimization. The 
comparison shows excellent agreement between the two 
programs. The temperature at which the equilibrium 
was calculated was 3000°K and the pressure was 10 
atmospheres . 

Methane Water Reaction . A comparison of 
results, with those of Oliver et al . (15), for the 
equilibrium reaction of five moles of Hydrogen and one 
mole of methane at 600°C and 1 atmosphere is 
presented in Table 4-9. Their results were obtained 
by the free energy minimization technique also. The 
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TABLE 4-8. Comparison of UDNH-RFNA Reaction Products at 
Equilibrium Predicted by Brandmair and 
Harnett and the Equilibrium Program at 
3000°K and 10 Atmospheres. 


Composition (Mole Percent) 

Species 

Brandmair and 
Harnett 

Equilibrium 

Program 

«2 

3.6569 

3.6567 

II 2 0 

40.7310 

40.7273 

OH 

9.5222 

9.5216 

CO 

7.3355 

7.3350 

co 2 

10.4574 

10.4563 

n 2 

23.2108 

23.2170 

NO 

0.9386 

0.9387 

N 

2.08 x 10' 3 

2.09 x 10" 3 

H 

.9554 

. 9553 

0 

.6136 

.6136 

0 2 

2.578 

2.578 
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TABLE 4-9. Equilibrium Composition for the Reaction 
of 5 H 2 O + CH 4 at600°C and one Atmosphere 
Pressure . 



Composition (Mole 

Percent) 

Species 

Oliver et al . 

Equilibrium 

Program 


CH 4 1.16 1.17 
H 2 43.36 43.35 
H 2 0 43.85 43.86 
CO 3.15 3.14 
C0 2 


8.48 


8.48 
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results show excellent agreement with chemical equilibrium 
composition calculated by the chemical equilibrium 
analysis program. 

Hydrazine Oxygen Reaction . The equilibrium 
composition of the . combustion of a stochiometric mixture 
of hydrazine and oxygen at 3500°K and 750 psia, as 
calculated by the equilibrium analysis program is present- 
ed in Table 4-10. These results are compared with those 
of White, et al ♦ (3) and Napthali (58, 59). White e t al . 
(3) showed that it took at least six iterations for the 
program to converge, our chemical equilibrium analysis 
program took three. The results with both authors 
shows excellent agreement. 

The previous five comparisons of the chemical 
equilibrium analysis computations with the results 
of othe*rs have shown the integrity of the program. A 
slightly modified version of the chemical equilibrium 
analysis program is used as a subroutine in the ABLATIN1 
and ABLATIN2 analysis programs to be discussed in 
Appendix A. The chemical equilibrium analysis program 
is discussed in detail in Appendix E. 

Free Energy Minimization: Is it the Best Method? 

Prior to 1958 all equilibrium calculations were 
performed using the equilibrium constant approach, the 
main proponent of which were Brinkley (9, 20) and Huff 
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TABLE 4-10. Equilibrium Composition of Products from 
the Combustion of a Stochiometric Mixture 
of Hydrazine and Oxygen at 3500°K and 
750 psia. 


Specie 

Starting 

Point 

White* 

Equilibrium 

Program 

Napthali** 

H 

.1 

0.040668 

0.0406415 

0.04118 

h 2 

.35 

0.147730 

0.147700 

0.14757 

h 2 o 

.50 

0.783153 

0.783222 

0.78358 

N 

.1 

0.001414 

0.01413 

0.001410 

N 2 

.35 

0.485247 

0.485249 

0.485210 

NO 

.1 

0.027399 

0.027395 

0.027480 

0 

.1 

0.017947 

0.017935 

0.017980 

o 2 

.1 

0.037314 

0.037303 

0.037570 

OH 

.1 

0.096872 

0.096848 

0.095820 


*Reference (3) 

**Reference (58, 59) 
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(17, 32, 23, 24). In that year, White, Johnson^, and 
Dantzig (3) of the RAND Corporation proposed to 
calculate chemical equilibrium of reacting gases by 
minimizing the free energy. Their method sparked 
considerable interest and controversy. It caused users 
to be divided into two camps: the faithful equilibrium 

constant formulators and the free energy minimi zers. 

The controversy flamed with claims and counter claims 
of assured convergence and of computational advantages 
offered by each method. As reported by Zeleznik and 
Gordon (11), "so heated became the controversy that 
when a panel discussion was arranged to discuss 
equilibrium computation in 1959, it was necessary to 
divide the panel into a free energy panel and an 

equilibrium constant panel." — 

In 1960 Zeleznik and Gordon (41) tried to settle 
the controversy by comparing the three general methods 
of Brinkley (9, 20), Huff (17 , 22 , 23 , 24) and Wliite 
with each other. They extended the form of these methods 
to permit condensed species as reaction products. An 
analytical proof was shown to indicate that the three 
methods could be placed in a form independent of the 
choice of components. Finally, they concluded that no 
one of the three methods offered "any significant 
computational advantages over the other two." However, 
in their most recent publications (11, 12), they seem 
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to reverse their standing by stating that free energy 
minimization does offer some computational advantages. 

Free energy minimization has some distinct 
computational advantages in that the number of 
equations to be solved is equal to the number of 
element plus the number of solid species plus one. 

For a carbon hydrogen system with one solid species, 
the number of simultaneous equations to be solved are 
four regardless of the number of gas species. This 
advantage cannot be claimed by the equilibrium 
constant formulation method. In addition there is 
no need to specify an independent set of chemical 
reactions as is required by the equilibrium constant 
formulation. This would indeed be quite a task in 

setting -up a problem for say a 100— species- and 2 

elements since a total of 98 independent chemical 
reactions must be postulated. For free energy 
minimization it is only necessary to estimate the 
number of species believed to be present at equilibrium 
and this greatly simplifies the laborious task of 
postulating a set of independent chemical reactions. 
Isomers are handled by free energy minimization as 
distinct chemical species while for the equilibrium 
constant formulation they require special treatment 
(59) . Initial estimates which are required to be 
rather close for equilibrium constant formulation (4) 
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are not for free energy minimization. It has been found 
from experience that when postulating a chemical system 
it is best to have all of the guessed input compositions 
within 2 orders of magnitude of each other. This 
reduces the number of iterations required for a specified 
temperature and pressure. 

All of the facts in this section tend to 
support the contention that as a general multipurpose 
scheme, for any number of species, free energy 
minimization is probably the best. This, of course, 
is not to say that existing programs implemented using 
the equilibrium constant approach are inferior. On 
the contrary, in the final analysis the success in 
implementing either technique is what determines its 

-usefulness. For -those- that would prefer to have a 

"feeling" for what reactions are affecting 
equilibrium the most, the equilibrium constant 
formulation may be used. This kind of information 
cannot be obtained directly from free energy minimiza- 
tion since it does not use chemical reactions in its 
implementation . 

When repetitive numbers of calculation are 
required, as in flow field calculation, free energy 
minimization might become time consuming and several 
specialized schemes (61, 62) have been developed to 
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take advantage of some special characteristic of the 
problem to reduce the time of solution. These schemes, 
of course, are quite restrictive for general 
applications; and they cannot be used for other than 
the solution to the particular problem for which they 
were developed. This is because any modification of 
the system, such as the inclusion of products not 
previously considered or the addition of a new element, 
would require a new analysis, algebraic manipulations 
and reprogramming. 

Recently, four additional techniques for 
computing chemical equilibrium composition for gases, 
have been published in the literature (63, 64, 65, 66) 
and they are briefly discussed. One is by White (63) 

■fn _ whrch— he— develops— a— technique- -b-a-sed— on— ' the— con--- 

tribution of each element to the total free energy of 
the molecule. In his new technique he makes claim 
that the method eliminates certain difficulties 
encountered by free energy minimization and in addition 
reduces the number of iterations required for con- 
vergence. The iterations are based not on the final 
value of the equilibrium composition of the chemical 
species, but on the final value of the Lagrange 
multipliers. There is, however, a drawback in that 
tables of the contribution, or partitioning of the free 
energy are necessary for this method, and these are not 
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yet available. Perhaps, as the author admits, not until 
these tables are readily available will there be much 
interest in this new method. The second technique 
is that of Lai, (64) which uses the mathematical duality 
between the minimum free energy and the mass conserva- 
tion equations to compute chemical equilibrium. Another 
dual technique is that of Passy and Wilde (65) and 
it is based on the duality between geometric programming 
and the minimization of the Gibbs free energy function. 
Several questions were raised during Lai's (64) 
presentation of his paper about the validity of his 
convergence scheme, and this raises the possibility 
of divergence in cases more difficult than the one 
presented by the author. 

-The— me thod— of— Ras sy_ and-Wi Ide— (-6.5-)— i s_ a 

completely new technique for solving single phase 
chemical equilibrium problems. The main advantage 
of transforming the problem to a posynomial 
minimization is a reduction in the number of variables 
to be solved. .Moreover, the number of. Lagrange 
multipliers is reduced to one which is interpreted 
physicially, as the total number of moles of the 
system at equilibrium. Since the method reduces the 
number of minimization variables, when implemented 
it will most probably reduce computing time. The 
method as derived is not applicable to polyphase 
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system and as such, its usefulness to this research 
where solid phases are considered is limited. The 
final and latest method in the literature is that of 
Meisnner et al . (66) . In contrast with the three 

previous methods which fall into the category of 
minimization schemes, the method of Meisnner and co- 
workers (66) is a variation of the equilibrium 
constant formulation which they call the "reactor 
series method". The method proposes that for 
calculating equilibrium compositions for a system 
in which many reactions occur, the system should be 
reduced to a series of individual reactions occuring 
separately. The main advantage claimed by the 
authors is the simplicity of the calculation which is 

predicated on— the fact that— since -only one— reaction 

is treated at a time simple algebra is involved in every 
step. However, this method suffers from the main dis- 
advantage of the equilibrium constant method, i.e., 
having to postulate a system of independent chemical 
reactions in comparison with free energy minimization 
where species believed to be present are the only things 
specified . 

None of the four methods reviewed, except 
perhaps that of Passy and Wilde (65) , which only 
applies to single phase system, offers any superior 
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advantage over the free energy minimization technique 

( 3 ). 

Summary 

This chapter has presented a detailed 
derivation of the free energy minimization technique. 
Comparison with other well known methods were made to 
illustrate that free energy minimization was more 
suitable for a number of reasons. The principle 
one being that the only thing required to compute 
chemical equilibrium is to specify the number of 
species rather than the number of independent 
reactions as is required for other methods. Four 
recent techniques which appeared in the literature 
were briefly discussed and none seem to offer any 
particular advantage over free energy minimization 
as it has been applied to this research. In the 
next chapter a detailed presentation of the non- 
equilibrium flow analysis is given. 
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CHAPTER V 


NON-EQUILIBRIUM FLOW ANALYSIS 

Introduction 

As mentioned in Chapter III, there are two limiting 
cases currently used to simplify the analysis of the flow 
of pyrolysis gases through the decomposition zone and char 
zone. These are to consider the flow frozen or in thermo- 
dynamic equilibrium. Frozen flow tends to under predict 
the energy absorbed in the decomposition and char zones 
since it does not account for the energy absorbed by the 
endothermic reactions . Equilibrium on the other hand , 
tends to over predict the total energy absorbed since it 
a ssumes that all the reactions proceed to their maximum 
extent. In order to obtain a more realistic analysis, a 

j 

third method was also studied in this research whereby the 
reaction rates of the chemical species are governed by 
finite reaction rates. This of necessity requires the 
determination of all important chemical reactions and the 
analysis of the appropriate kinetic data corresponding to 
these reactions. This method of analysis, therefore, 
requires the laborious task of collecting, cataloguing and 
evaluating reaction kinetic data as it appears in the 
literature . 
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A procedure was developed to help simplify, but not 
totally eliminate, the difficulty encountered in evaluating 
kinetic data for the reactions, and to help narrow down the 
choice of selecting information from sources which some- 
times do not agree among themselves. This procedure is 
called the isothermal analysis- technique, and it is ex- 
plained in a latter part of the chapter. In addition to 
the painstaking data handling and analyzing procedure there 
is the added task of numerically solving the equations. 

The non-equilibrium flow analysis is an order of magnitude 
more complex than the chemical equilibrium analysis of 
Chapter IV. The latter involves the solution of a set of 
algebraic equations with the energy equation; while the 
non-equilibrium flow analysis requires the solution of a 
set of non-linear ordinary differential equations (species 
continuity) with the energy equation. The numerical solu- 
tion of these equations can become stiff (10) as will be 
explained subsequently. 

Mathematical Formulation 

The phenomenological equations that permit the calcula- 
tion of the reaction rates for an arbitrary number of 
simultaneous chemical reactions involving an arbitrary number 
of chemical species is presented in this section. 
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A chemical reaction can be written in general as: 




A 

j 


i 


1 . 


. .L 


(5-1) 


For this ith chemical reaction the r . . and p . . represent 
the stochiometric coefficient of the reactants and products 
respectively. Aj represent the chemical species j. The 
forward and reverse reaction rate constants are k*. and 

r i 

k r ^, respectively. Assuming that there are a total of 1 
chemically reacting species whose composition can be 
determined by L chemical reactions, the rate of reaction 
of the j species , Rj , can be described by the following 
phenomenological equation : 



I 

1 r . . 

(Pi-j-r.J'tk H C 13 
13 ^ f i k=l K 



1 

n 

k=l 



(5-2) 


where C R is the concentration of component K, and k r ^ 

are the forward and reverse reaction rate constants, and 

r^j and pij are the power on the concentration of the 

> » 

reactants and products respectively, rij and pij are 
equal to their corresponding stochiometric coefficients for 
elementary reactions only. 

Equation (5-2) is a convenient and general formulation 
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for expressing the reaction rate of the jth species in 
L simultaneous chemical reactions. This form of the rate 
equation for species i is readily used in the computer 
implemented numerical solutions of the non-equilibrium 
flow analysis. The stochiometric coefficients stored in 
matrix form, where the species are identified by their 
row location and the reaction by their column location. 
Since there is no algorithim to relate the rows and columns 
for use in a generalized, all purpose input subroutine, 
i. e. , the matrix of coefficients is a sparce matrix. The 
forward and reverse reaction rate constants are convenient- 
ly expressed in the following functional form: 


k i = A i T " Si exp <- E i/ RT ) i=l-..L (5-3) 

when S i =»0 Equation (5—3) reduces to the well known 
Arrhenius expression. Ai is the frequency factor and 
is the activation energy. 

Equations (5-2) and (5-3) are simple to program. 
However, the selection and analysis of the data for the 
important chemical reactions and the subsequent numerical 
problems that are encountered in describing chemically 
reacting flow, makes this analysis more difficult than 
the equilibrium case. The selection of the chemical 
reactions and the numerical problems encountered in their 
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solutions are explained in the next two subsequent sec- 
tions . 

Criteria for Reaction Selection 

When developing a kinetics analysis for a chemically 
reacting flow system one must first determine what are 
all the possible chemical reactions. That is, one must 
perform an exhaustive analysis of all possible chemical 
reaction combinations and select from these, as a first 
approximation, those reactions which are thermodynamically 
feasible. To determine all the possible chemical reactions, 
three main sources were used (2, 3, 4). The work by 
Hockstein (2) covers over 20,000 possible reactions, not 
all applicable to this research. However , it has ari“ex- 
tensive list of reactions for the carbon-hydrogen and 
carbon-oxygen systems. Pike (4) also covered most of the 
important hydrocarbon reactions which were of particular 
interest to this research. Moreover, Bahn (3) recently 
listed all of the possible reactions in the H-O-N system. 

The first step in selecting the appropriate reactions 
was to eliminate those reactions which had reactants and 
products which could not be formed by any combination 
from the species present in the pyrolysis gases. A de- 
tailed explanation on how these compositions were obtained 
is given in Appendix F. 
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The next step was to determine those reactions from 
the remaining set which were thermodynamically feasible. 
Those that were feasible, were further analyzed using the 
isothermal analysis technique which is explained later in 
this section. 

r 

Composition of Pyrolysis Products : It has been 

mentioned that the composition of the pyrolysis gases has 
an important effect in the reaction selection process. 

It eliminates from consideration those reactions whose 
reactaBts could not be formed by any combination of the 
pyrolysis products. Thus, it is most important to estimate 
the com position of the p yrol ys is gases with as much cer- 
tainty as possible. Note that the composition of pyro- 
lysis gases only affects the non-equilibrium analysis since 
equilibrium only requires the elemental composition to be 
known, and these are known from the chemical formulae of 
the components of the ablative composite and their composi- 
tions . 

First attempts to study the non-equilibrium flow of 
pyrolysis products relied on the composition obtained from 
two sources. One being the composition obtained assuming 
chemical equilibrium (5) and the other by the analysis of 
the degradation products of low density nylon phenolic 
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resin by pyrolysis gas chromatography (6, 7) . As noted 
in Appendix F, the unavailability of accurate analytical 
procedures and thermophysical properties for the high 
molecular weight pyrolysis products ( i. e., phenol, cresol, 
toluene, etc.) left a region of definite uncertainty as to 
the composition of pyrolysis products. As a result, the ma- 
jor components of the degradation of the pyrolysis products 
were identified as methane, hydrogen, carbon dioxide, carbon 
monoxide and nitrogen by many (8, 9), with unknown quan- 
tities of water and high molecular weight residue completing 
the analysis. 

Subsequent research by Sykes (10) confirmed the presence 
of phenol-based materials as primary constituents in the 
high - molecularweight residues. — Table— 5—1— lists-the-species— 
and the estimate of the composition of the pyrolysis products 
resulting from the degradation of a 40 percent (by weight) 
nylon, 60 percent phenolic resin ablative composite. How 
these estimates were arrived at is presented in Appendix 
F. 

Equilibrium Conversion of Reactions : The second step 

in the development of a realistic kinetic analysis is to 
determine from the list of possible reactions which are 
thermodynamically unfavorable. By computing the equilibrium 
conversion of the reactions over the temperature range of 
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TABLE 5-1 Estimate of the Representative Composition 
of the Pyrolysis Products for a 40 percent 
Nylon, 60 percent by Weight Phenolic Resin 
Ablative Composite. 


Species 

Mass Percent 

Mole Percent 

H 2 

2.57 

20.92 

ch 4 

3.86 

3.90 

c 2 h 2 

3.88 

2.41 

c 2 h 4 

3.88 

2.24 

c 2 h 6 

0.64 

.35 

C 6« 6 

2.57 

.53 

CgHgOH 

23.10 

3.97 

CO 

4.65 

2.40 

CO 2 

4.60 

1.69 

H 2 o 

7.29 

6.45 

n 2 

3.95 

2.20 

C (solid) 

39.01 

52.94 



Total 


100.00 
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interest (500°F - 5500°F) one could examine the beha- 
vior of the equilibrium conversion of each reaction with 
temperature. The rule of thumb used was to eliminate 
from further consideration those reactions which had an 
equilibrium conversion of 5% or less at 5500°F. This is 
justifiable since kinetically their conversion would 
have been much less. We are assuming, of course, that 
the equilibrium conversion is determined in the endother- 
mic direction. In Appendix H, the equilibrium conversion 
with respect to temperature for 99 possible chemical 
reactions is presented. These reactions were obtained 
from possible combinations of the pyrolysis products as 
presented in Table 5-1. It should be noted that many of 
the reactions considered, a lthough ther^dyhamica lly 
feasible, were eliminated, since either one of the reacting 
species did not exist as a pyrolysis product, or could not 
be formed by any combination of important reactions. One 
example of this is the dissociation of N 2 0 as given by 
reaction 7 in Table H-3. As shown in Table 5-1, there is 
no N 2 0 present as a product of pyrolysis degradation, 
and practically none can be formed by the oxidation of 
nitrogen as indicated by the equilibrium conversion of 
reaction 17 of Table H-3. Therefore, there was no need to 
consider the dissociation of N 2 0 in the kinetic analysis, 
and this sort of reasoning was used to eliminate other 
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reactions. 

Isothermal Analysis Technique ; In Tables H-l through 
H-7 a list of possible chemical reactions is presented 
along with the equilibrium conversion of these reactions. 

In these tables it is shown that many of the possible 
reactions do not have appreciable conversions , even at 
the higher temperatures and could be, therefore, elimina- 
ted from further consideration. In many instances as 
previously mentioned, reactions which were thermodynami- 
cally favorable were eliminated because their reactants 
did not exist, or could not be formed as products from 
the pyrolysis reactions. This point was illustrated with 
the N 2 O dissociation reaction. However, there was a 
remaining group of reactions (too many to have been 
considered in the analysis) that could not be eliminated 
from consideration by any of the two previously mentioned 
procedures. Therefore, a criteria was developed to determine 
the kinetic importance of the remaining reactions in the 
temperature range from approximately 500°F to 5500°F. 

By analyzing the conversion of the flow of an equal molal 
mixture of the reactants, one was able to infer the im- 
portance of such a reaction. The mass flux used was 
0.01 lbs/ft^-sec. and a char thickness of 0.25 inches 
with a porosity of 0.8. 

In Figure 5-1 the conversions are presented for two 
reactions. These are the thermal decomposition of ethylene 



Percent Conversion 


Conditions 

Mass flux =0.01 lb/ft -s< 
Pressure=1.0 atm / 
Porosity=0.8 / 
Thickness=0.25 in. / 


Ethvlene 


Acetylene 


C 2 H 2 =2C + H 2 


500 1000 1500 2000 2500 300 
Temperature °F 

Figure 5-1. Conversion of Pure Ethylene and 
Pure Acetylene at Isothermal Con 
ditions (500°F - 3000 F). 
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and acetylene species. It can be concluded by examining 
this figure that there could be no conversion of ethylene 
within the char for temperatures less than 1000°P and 
no conversion of acetylene for temperatures less than 
2000°F. However, the thermal cracking of ethylene would 
be of importance for temperatures above 1000°F and for 
acetylene above 2000°F. In most instances the other 
reactions examined remained frozen below 1200 to 1500°F. 

But usually once the reaction started it would react very 
fast. This behavior, which is typical of reactions with 
high activation energies {>50 Kcal/mole) and high frequency 
factors (>10 15 ) , was observed in the analysis. This very 
fast take-off of many of the reactions considered in the 
kinetic analysis was the cause of most of the numerical 
problems encountered in the integration of the equations 
of change. This behavior manifested itself in a condition 
known as stiffness, and more is said about this in a 
latter section. 

To make the discussion quantitative, the conversion 
of a reactant is defined as the ratio of the amount 
consumed by reaction to the amount initially present. To 
determine the conversion of a chemical reaction a material 
balance is made on component j flowing through a volume- 
tric section of the char having a cross-sectional area, 
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A, and a width A 2 as shown in Figure 5-2. If N. is 

3 

the molal flux of component j at z, the material 

balance on component j for steady flow in the 2 direction 

is: 


(N.M A) + (R M .A) . Az - (N.M .A) A = 0 
3 wj z ] W] j w] ' z+Az 


(5-4) 


dividing the above by M^A results in the following 
equation: 


(N_.) 


z+Az 


- v* 


(5-5) 


The above equation expresses the change in molal 
flux of each species, as it reacts through th e is othermal 
reactor. This equation is solved by a simple finite 
difference scheme from which the conversion Xj is cal- 
culated as follows: 


X. = (N. 

3 30 


- N.)/N. 

3 30 


(5-6) 


where Nj Q is the molal flux of component j entering the 
char zone. By using this scheme the set of important 
chemical reactions was obtained and they are listed in 
Table 5-2. The computer program implemented solution 
for the finite difference Equation is listed in Appendix 
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TABLE 5-2 Chemical Reactions of Importance Considered 
in the Non-Equilibrium Flow Analysis 


1. ch 4 = ch 2 + h 2 

2. CH 4 = CH 3 + H 

3. 2CH 3 = C 2 H 6 

4. C 2 H 6 = C 2 H 4 + H 2 

5. C 2 H 4 = C 2 H 2 + H 2 

6. C 2 H 2 = 2C + H 2 

7. C 2 H 2 = C 2 H + H 

8. C + H 2 0 * CO + H 2 
(s) 

9. C + CO, = 2C0 

(S) “ 

10. H 2 + M = 2H + M 

11. H 2 0 +M=H+OH+M 

12. H + C02 - CO + OH 

13. C0 2 = CO + 0 

14. CgHgO + H 2 = H 2 0 + CgHg 

15. CgHg = 3C 2 H 2 
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D. The chemical reactions selected to simulate the 
chemical behavior of the pyrolysis gases are presented 
in Table 5-2. 


Selecting the Right Data 

The procedure developed to analyze and screen the 
many possible chemical reactions which occur in the 
C-H-O-N as explained in the previous section, was helpful 
in selecting a set of reactions which would most closely 
represent the actual chemical behavior of the pyrolysis 
gases. We recall that this selection is based on the 
use of the isothermal analysis technique. This technique 
computed the isothermal rate of conversion of a reaction 
over^a“temperature — ranges — The degree of conversion over 


this range allowed us to form a judgement as to the im- 
portance of a given reaction. If a reaction was judged 
to be important (conversion> 5 percent) it was included in 
the analysis, otherwise it was eliminated. Once the reac- 
tions were selected, one additional problem remained. 

This was the selection of data from a number of sources 
which qualitatively predicted the right trend in the reac- 
tion but which predicted quite different conversions. In 
other words, the problem was what to do in cases where two 
sources of data existed for a given reaction, which predict 
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different conversions at given isothermal conditions. The 
approach to this problem is explained below. 

Comparison of Isothermal Conversion With Equilibrium 
Conversion 

The comparison of the isothermal conversion with the 
equilibrium conversion over the temperature range of 500°F 
to 5500°F eliminated a great deal of the data selection 
problem. In most cases, the comparison showed that only 
one of the data sources would remain within the equilibrium 
constraints for a given reaction. Therefore, the choice 
of which source of data to use was obvious in these cases. 
However, there were instances in which the reaction rate 
data predicted concentrations in violation of known ther- 
modynamic equilibrium constraints. This is attributed to 
the fact that in most cases the kinetic data was being 
extrapolated beyond the temperature range of applicability. 

Rather than to eliminate these reactions it was 
decided to force them to be within the equilibrium cons- 
traints. This was achieved by using a reverse reaction 
rate constant calculated from the well known equation which 
relates the reverse constant with the equilibrium constant 
and the forward reaction constant. In doing so, a set of 
reactions that met equilibrium constraints was obatined. 
Again more than one source of data existed for each one of 


these reactions 
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Fortunately, of the 15 reactions that were selected to 
simulate the kinetics, only four fell in this category, 
and are shown in Table 5-3. Therefore, a parametric stu- 
dy was conducted to evaluate the effect that each of 
these four reactions had on the kinetics analysis when 
a different source of data was used. The results showed 
less than one percent variations between the predicted 
values of the energy absorbed for each reaction. These 
small variations are attributed partly to the use of 
a reverse reaction rate constant (consistent with equi- 
librium) which tended to dampen out the differences in 
the values of the forward reaction rate constants. An 
explanation on how the parametric studies were conducted 
is given in the following chapter. 

Discussion of the Kinetic Data for the Important Chemical 
Reactions in the Char Zone 

Table 5-4 presents the list of the 15 reactions used 
in the non-equilibrium flow analysis. It also presents 
the kinetic data and literature sources of these reactions 
The source of the preferred value is underlined. 

Methane Decomposition : Reactions 5-1 and 5-2 were 

used to describe the decomposition of methane. Reaction 
5-1 decomposes to CH 2 and H 2 . 


Data for this reaction were 






TABLE 5-3 Reactions for which more than one set of 
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reported by Palmer and Hirt (36) and Kozlov and Knorre 
(37) . An activation energy of 91 Kcal/mole and a fre- 
quency factor of 4.5 x 10 ^ sec“^ were reported by Palmer 
and Hirt (36) . Kozlov and Knorre (37) reported an acti- 
vation energy of 101 Kcal/mole and a frequency factor of 
1.26 x 10^ 4 sec"*. The data from both of these sources 
were used in the prediction of the total energy absorbed. 

It was observed that there was a negligible difference 
(less than one half of one percent ) between the predicted 
values of the energy absorbed using the data of Palmer 
and Hirt (36) and that of Kozlov and Knorre (37) . As a 
matter of interest to the reader , the calculated values 
reported in the latter part of this chapter for the total 

energy absorbed were those us ing thedata o f — Ko zlovand 

Knorre (37) . The reverse reaction rate constants were 
not reported in the literature, however, they were easily 
computed knowing the forward reaction rate constant and 
the equilibrium constant. The equilibrium constant was 
calculated as a function of temperature using a modified 
version of the equilibrium analysis program described in 

Chapter IV. These values were fitted using an expression 
( A+B/T) 

of the form e' The values for the constants A and 

B are given in Table 5-5 for the 15 chemical reactions used 
in the kinetic model. 

A second competing reaction was used to describe the 
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decomposition of methane. This is reaction 5-3 which 
shows the decomposition of methane to methyl radical 
(CH 3 ) and hydrogen radical (H) . For this reaction the 
data published by Palmer and Hirt (38) was used. The 
activation energy reported was 103 Kcal/mole with a fre- 
quency factor of 1.0 x 10^ sec - ^. This data was reported 
by Pike (39) as preferred in his critical evaluation of 
rate data for reactions in the C-H-O-N system. The range 
of activation energies in the literature varied from a 
low value of 73 Kcal/mole to a high of 103 Kcal/mole (36, 
38,40,41,42,43). Frequency factor variations between 
4.5 x 1013 to 1.0 x lO ^ 5 sec -1 were also reported. 

Formation of Ethane : React i on 5-3 wa s used to re- 
present the formation of ethane from two methyl radicals. 
These radicals are formed in the decomposition of methane. 
The activation energy used was 28.2 Real/ (gm. mole-sec) . 
There was only one source of data known to this author 
for this reaction and it is reported by Steacie (44) . 

Decomposition of Ethane : Reaction 5-4 was used to 

describe the decomposition of ethane. In this reaction 
ethane decomposes to ethylene and hydrogen. Data for this 
reaction have been reported by a number of authors. These 
include Gulyaev and Polack (45) , Kozlov and Knorre, Shah 
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(46) , Bartlett and Bliss (47) and Steacie and Shane (48) . 

An activation energy of 69 Kcal/gm.mole and a frequency 
factor of 9x10^3 sec - ^- were reported by Gulyaev and Po- 
lack (45) . Kozlov and Knorre (41) reported an activation 
energy also of 69.0 and a frequency factor of lxlO 14 sec”*. 
This is very close to the data by Gulyaev and Polack (45) 
which reported a frequency factor of 0.9xl0^ 4 sec"^ and 
an activation energy of 69.0 Kcal/gm.mole. An activation 
energy of 64.1 Kcal/gm.mole and a frequency factor of 
3.14x10^ sec”^ were reported by Bartlett and Bliss (47) 
which is also close to the data of Kozlov et.al. (41) and 
Gulyaev et.al. (45). Shah (46) reported an activation 
energy of 83 Kcal/gm.mole and a frequency factor of 6.04x 
1016 sec~l. Both the acti vation energy and the frequency 
factor of Shah are way out of line to those of the three 
previously mentioned sources. For this reason the data 
of Shah was not used. 

The range of activation energies and frequency factors 
reported (41, 45, 46,47,48) were 64 to 82 Kcal/gm.mole and 
3.14xl0 13 to 6.04xl0 16 respectively. 

The calculations of the total energy absorbed reported 
in the latter part of this chapter were based on the data 
of Kozlov and Knorre (41) for ethane decomposition. 
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Decomposition of Ethylene : The thermal degradation of 

ethylene to acetylene and hydrogen is given by reaction 
(5-5) . Data for this reaction has been reported by Kozlov 
and Knorre (41) , Gulyaev and Polack (45) and Shah (46) . 

The activation energy and frequency factors reported by 
Kozlov and Knorre (41) and Gulyaev and Polack (45) were 
identical. Both reported an activation energy of 40.0 
Kcal/gm.mole. However, Kozlov and Knorre (41) gave a tem- 
perature range of 1300-2000°K where the data would be 
applicable, while Gulyaev and Polack (45) gave a temperature 
range of 1600-3700°K. The range of activation energies 
and frequency factors reported (41,45,46) were 40.0 to 
76.0 Kcal/gm.mole and 2.57x10® to 1.8xl0* 3 sec 1 , respec- 
tively. 

A value of 40. Kcal/gm.mole and a frequency factor of 
2.57x10® sec“"l were used in the kinetic analysis. 

Acetylene Decomposition : The thermal decomposition 

of acetylene to carbon and hydrogen is given by reaction 
(5-6) . The number of data sources for this reaction was 
greater than average. One would think that a reaction so 
well studied, and of such commercial importance would show 
good agreement among the sources. Unfortunately, the 
opposite was the case. The greatest discrepancy in the cal- 
culated values of the energy absorbed were observed when 
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using different sources of data for this reaction. The 
difference amounted to about four percent. 

Data for this reaction were reported by Kozlov and 
Knorre (41) , Gulyaev and Polack (45) , Shah (46) , Happel 
and Kramer (49), Aten and Greene (50), Leroux and Mathieu 
(51), and Chase and Weingberg (52). Both Kozlov et. al. 

(41) and Gulyaev et„al. (45) reported a frequency factor 
of 1.7x10® and an activation energy of 30.0 Kcal/gm.mole. 
Kozlov et. al. (41) gave 2000°K as the temperature at which 
these data were applicable. Gulyaev et. al. gave a much 
wider temperature range of 1600°K to 3700°K. 

Shah (46) reports an activation energy of 62 Real/ 
gm.mole, twice that of the two previous authors, and a 

frequency_f actor of 9.7xl0 1 ^ which is50 ,000 times greater. 

Because of the higher activation energy and larger frequency 
factor, the data of Shah (46) would predict acetylene de- 
composition over a smaller range than the data of Kozlov 
(41) and Gulyaev (45) . The data of Shah for ethane and 
ethylene decomposition also poses these characteristics. 

The frequency factor of Happel and Kramer (49) , 
S.lxlO 10 , is of the same order of magnitude as that repor- 
ted by Shah (46) . However, the activation energy is 
approximately three and a half times smaller. 

Aten and Greene (50) report an activation energy of 
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of 14.5 Kcal/gm.mole, slightly lower than Happel and Kramer 
(49) , but half that of Kozlov (41) and Gulyaev (45) . A 
frequency factor of 4.0x10^ was given, which is twenty- 
three times greater than Kozlov and Gulyaev but 3000 
times smaller than the one reported by Happel and Kramer 
(49) . It is also about 2500 times smaller than the data 
of Shah (46) . The temperature range for which the data 
is recommended is 900-1700°K. A much smaller temperature 
level than the 1600-3700°K range recommended by Gulyaev 
(45) . 

Leroux and Mathieu (51) report a frequency factor of 
1.0x10^ and an activation energy of 12.5 Kcal/gm.mole. 

These values are completely out of line with those repor- 
ted in the literature by-other authors.- The temperature 
range for which these values are recommended are 298-398°K. 
Obviously, too low for the temperature range of interest in 
this research. The values of Leroux and Mathieu were never 
considered and they were not used in the study of the effects 
of kinetic data on the predicted values of the energy 
absorbed. The basis for this decision is obvious. 

The data of Kozlov and Knorre (41) and Gulyaev and 
Polack (45) were used for reaction (5-6) . 

Reaction (5-7) was also used to describe the thermal 
decomposition of acetylene. The data of Eschenroeder and 
Lordi (54) was used to simulate this reaction. They reported 
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an activation energy of 35.5 Kcal/gm.mole and a frequency 
factor of 4.5X10 11 . This was the only single reference 
found for this reaction. 

Carbon-Steam Reaction : The reaction of carbon with 

steam to produce carbon monoxide and hydrogen is given by 
Equation (5-8) . This reaction has been extensively 
studied and the amount of data available is copious (55, 
56,57,58,59,60,61). The reaction has been established to 
be first order and activation energies of 26-90 Kcal/gm.mole 
have been reported. Most of this data in the literature is 
geared towards industrial applications and therefore the 
form in which the data was presented was not easily ap- 
plicable to a generalized reaction kinetics routine. 
Fortunately, the data of Walker et.al. (61) for the graphite 
steam reaction was presented in Arrhenius forms, which 
conforms to the type of expression used in the kinetics 
subroutine. Walker et.al. (61) report an activation energy 
of 71.9 Kcal/gm.mole, and a frequency factor of 2.3X10 11 
sec -1 . The activation energy is close to the value of 70 
Kcal/gm.mole reported by Pike (4) in his literature 
evaluation. 

Carbon -Carbon Dioxide Reaction: This reaction has also 
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been thoroughly studied in the literature because of 
important industrial applications. Data for this reaction 
has been reported by Wu (64) , Gadsby et. aU (65), Lewis 
and co-workers (66) , Harper (67) , Austin and Walker (68) , 
and Glovina (69) . The form of the reaction rate for the 
data of references (28-31) is: 

r = k l P co 2 

co 2 (5-7) 

+ k 2 P co + k 3 P co 2 ) 

where k^, k 2 and k 3 are constant and P C q 2 and P C0 are tlie 
partial pressure of C0 2 and CO respectively. The previous 
form of the equation was not used. Instead the form used 

was~k]_ (A)— where— (A)— represents— the-concentration-of—G0 2 -i 

In other words, the reaction was assumed to be first order 
with the concentration of CC> 2 . This approach was taken 
due to a suggestion by Swann (70) in which he pointed out 
that the C0 2 data available in the literature was for partial 
pressures of C0 2 much lower than the ones that would be 
encountered in re-entry simulation. Therefore, the litera- 
ture data would not be any more accurate than the first 
order approximation. This approach was used by April (71) 
and Pike et. al. (72). They reported an activation energy 
of 50 Kcal/gm.mole, and a frequency factor of l.OxlO 6 . 
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However, the activation energy used in our research was 
35 Kcal/gm.mole and a frequency factor of 1.0x10*^. This 
is the data reported by Walker (62) . 

Hydrogen Dissociation : The hydrogen dissociation 

reaction. Equation (5-10) , is another reaction that has 
been thoroughly studied in the literature. Ellis et. al. 

(73) , Kaskan and co-workers (74) , Fenimore (75) , and 
Gardiner and Kistiakowsky (76) among others, have reported 
values for the activation energy and frequency factor of 
hydrogen dissociation reaction. The data of Kaskan and co- 
workers (74) was used in the kinetic analysis. They re- 
ported an activation energy of 103.2 Kcal/gm.mole and 

18 

a frequency factor of 3.6x10 . Similar values were also 

reported by Bartlett and Bliss (45) in their studies of 
methane decomposition. 

It was observed in the kinetic analysis that no 
appreciable dissociation of hydrogen occurs below 2500°F 
which is in line with the high activation energy reported 
in the literature. This is also consistent with thermo- 
dynamic equilibrium calculations obtained with the chemical 
equilibrium program. This reaction was not used by April 
et.al. (71, 72) since their analysis was restricted to a 
maximum temperature of 3000°F, and is in line with our 
findings that reactions do not appreciably occur below 
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temperatures of 2500°F. 

Water Dissociation : The thermal dissociation of water 

to hydrogen and hydroxyl radical is given by reaction (5-11) . 
Harper (67), Ellis et.al. (73) and Kaskan (74) report values 
for the activation energy and frequency factor of water 
dissociation. Harper reports an activation energy of 7.7 
Kcal/gm.mole and a frequency factor of 2 ,5xl0^ 2 cm^/mole-sec. 
This is in marked contrast with higher values reported by 
Ellis et.al. (73) and Kaskan (74). The values used in the 
kinetic analysis were those reported by Ellis et.al. (73) 
and are 103.1 and 3.6xl0 x for the activation energy and 
the frequency factor respectively. 


Carbon Dioxide Dissociation ; The carbon dioxide dis- 
sociation is given by reaction (5-12) . Data for this reac- 
tion has been reported by Harper (67), Jensen and Kurzuis 
(78), Mahan and Solo (79), Davies (77), and in the literature 
review of Kaskan and Browne (74) . Pike (4) in his literature 
review recommends the values reported by Davies. These are 
an activation energy of 74.4 and a frequency factor of 
2.44X10 11 . The values used in the kinetic analysis were 
reported by Harper (67) and are: 71.9 Kcal/gm.mole for the 
activation energy and 2.3x10*^ for the frequency factor. 
Uncertainties in the value of the kinetic data for carbon 
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dioxide dissociation could not have had a very large impact 
in the predicted value of the energy absorbed since the ini- 

_ 9 

tial mole fraction of CO 2 is of the order of 10 . 

Carbon Dioxide-Hydrogen Reaction : The carbon dioxide 

hydrogen reaction which results in the production of carbon 
monoxide and hydroxyl radicals, is given by reaction 5-13. 
This reaction is important only at higher temperatures 
(>3000°F) where the concentration of the active species 
(hydrogen and OH radicals) are significant. An activation 
energy of 33 Kcal/gm.mole and a frequency factor of 1.3 x 
10 15 reported by Fristrom (78) were used in the kinetic 
analysis. This reaction was added to the kinetic analysis 
for the sake of completeness. Its contribution to other 
energy absorbing reactions is insignificant. 


Phenol Hydrogenation ; The phenol hydrogenation reaction 
to water and benzene is given by reaction (5-14) . This 
reaction was included in the analysis based on the experi- 
mental studies of April (71,72). Phenol reactions have been 
studied by many investigators (81, 82, 83, 84, 85). However, 
no kinetic data has been reported for this reaction. April 
(71) calculated a frequency factor of 2xl0 13 cm 3 mole -1 sec -1 
based on kinetic theory. He recommended an activation ener- 
gy of 45 Kcal/gm.mole based on his comparisons with heats of 
reactions of phenol related components, such as ethyl benzene 
(48.9 Kcal/gm.mole). O-Xylene (47.3 Kcal/gm.mole), mesitylene 



226 


(47.6 Kcal/gm.mole) and hydrindol (45.8 Kcal/gm.mole) . 

Benzene Decomposition : The decomposition of benzene 

is given by reaction (5-15) . The data for this reaction was 
reported by Pike (4) . An activation energy of 52 Kcal/gm- 
mole and a frequency factor of 1.4x10^ sec""* was used. 

Decomposition Kinetics of Ablative Composites 

In the earlier part of this chapter we discussed the 
kinetic of the pyrolysis gases in the char zone. In this 
section we will briefly explain the technique and the data 
used to describe the kinetics of the degradation of the po- 
lymer composite. 

The analysis of decomposition in depth requires that the 
rate of mass loss, or rate of temperature change with densi- 
ty be known. Data for the densities of virgin and degrading 
composites as a function of temperature have been reported 
by Sykes and Nelson (86) and Madorski (82) . 

The data of Sykes and Nelson (86) is particularly useful 
since it is for phenolic-nylon resins; The data was obtained 
using thermogravimetric analysis techniques. In Chapter II, 
the importance of using a kinetic equation of the Arrhenius 
type to correlate the experimental mass loss rate data for 
polymers was discussed. It was indicated in this chapter 
that the mass loss of material was affected by the heating 
rate, and that this effect was a source of difficulty for 
earlier researchers modeling the decomposition process. 
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Sykes and Nelson (86) used a pseudo-order kinetic ex- 
pression of the Arrhenius type to eliminate this influence. 
The equation is: 


dp^ n . 

= p i,o {(p i" p c ,i )/p i,o } 1 A i exp(-E./RT) 


(5-8) 


dt 


Equation (5-8) expresses the rate of change of density 
of a polymer with temperature, p. is the initial or virgin 
density. (* is the density of the material at temperature 
T, p c the residual or char density, and A^ and Ej_ are the 
well known frequency factor and activation energy parameters. 

In the present mathematical analysis it has been assumed 
that when an ablative composite degrades, it degrades inde- 
pendently of the other components ; that is , no interaction 


is assumed to occur among the composites. There is no known 

experimental data that takes into account these interactive 
effects. Therefore, of necessity, this simplifying assump- 
tion has been made. This is mathematically expressed as: 


dp _ 
dt 


v dP 
dz 


(5-9) 


and for quasi-steady flow the time dependent term can be 


modified by: 


dp 

dt 


v dp 
dz 


(5-10) 
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This , knowing the kinetic parameters for the various com- 
ponents in a blend of virgin material, specifying the surface 
recession velocity, and the temperature history from the 
energy equation, the variation in density of the virgin ma- 
terials can be predicted by the use of Equations (5-8) , 

(5-9) and (5-10) . 

Differential Thermal Analysis : The DTA thermograms of 

nylon, phenolic and phenolic microballoons which compare the 
phenolic-nylon composite as obtained by Sykes and Nelson 
are shown in Figure 5-3. The results shown were obtained by 
heating the material at 10°C/min in a helium atmosphere to 
the temperature shown. The reactions which occurred upon 
heating are either endothermal or exothermal, with the endo- 
thermal reactions extending downward from the base line 
A T=0^ ~ 

As shown in Figure 5-3 the thermogram of nylon undergoes 
two endothermic processes, one with the peak at 260°C and 
the second at approximately 417°C. As explained by Sykes 
and Nelson (86) , the first endotherm is associated with the 
melting of the Crystalline portion of the polymer (heat of 
fussion = 73.6 KJ/Kg of original material) while the second 
corresponds to decomposition. The decomposition, which 
occurs between 350°C and 500°C, absorbs 630 KJ/Kg of original 
material. According to Sykes and Nelson (86) "the apparent 
exothermal portion of the thermogram after 470°C results from 
the change in the total heat capacity accompanying decomposi- 
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TEMPERATURE, ° C 


Figure 5-3. Differential Thermal Analysis Thermogram 

of Nylon, Phenolic and Phenolic Micro- 
balloons as Reported by Sykes and Nelson 
( 86 ). 
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tion”. 

The center thermogram of Figure 5-3 is for phenolic 
which is used as the binder in the phenolic -nylon composite. 
The thermogram shows an exothermal reaction at 265°C fo- 
llowed by two overlapping endothermal decomposition reactions 
at 450°C and 625°C. The exothermal reaction evolves 68 KJ/ 
Kg, and the pyrolysis reactions between 350°C and 850°C 
absorbed 293 KJ/Kg. 

The thermogram for the phenolic microballoons, which is 
shown in the lower portion of Figure 5-3, is composed of 
small hollow micro-spheres of phenolic. The decomposition 
of this material is very similar to the phenolic binder. 

The exothermal reaction occurs approximately at 270°C and 
evolves about 47 KJ/Kg while the two overlapping endothermal 
reactions between 350°C and 850°C absorb 377 KJ/Kg. 


Thermogravimetric Analysis ; In conjunction with their 
DTA work, Sykes and Nelson (86) report the results of the 
thermogravimetric analysis (TGA) of the three separate cons- 
tituents. Figure 5-4 shows the TGA plot of the residual mass 
fraction as a function of temperature. 

The TGA data of Sykes and Nelson (86) shows one sharp 
mass loss between 350°C and 500°C for nylon. Almost all of 
the nylon is converted to gaseous product through this tem- 
perature range with less than ten percent remaining at 500°C. 

The TGA thermogram of phenolic and phenolic micro- 
balloons also shown in Figure 5-4, shows the mass loss rate 
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TEMPERATURE, °C 


Figure 5-4. Thermogravimetric Analysis Thermogram 

of Nylon, Phenolic and Phenolic Micro- 
balloons as Reported by Sykes and 
Nelson (86). 
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to continue through 850°C. Sykes and Nelson (86) determined 
that the residue of both phenolic materials was 54 percent 
of their original mass. 

The results of analysis of the three materials used in 
the composite phenolic-nylon material were very conveniently 
summarized by Sykes and Nelson (86) in a form which could be 
conveniently used by Equation (5-8) . Table 5-6 shows their 
summary. 

The data shown in Table 5-6 was used in this research 
to compute the mass loss rate of the phenolic-nylon resin. 

Silicone Elastomers : The decomposition of silicone 

elastomers was computed using the same procedure as for 
phenolic-nylon: i.e., using Equation (5-8) to compute the 

degradation process. 

At the time these computations were performed, Sykes 
had not published his results and the numbers presented in 
Table 5-7 were obtained by private communication with him 
( 88 ) . 

As shown in Table 5-7, a single reaction was used to 
describe the decomposition of this elastomer. 

The advantage of using Equation (5-8) as a general 
equation is that it will allow future use of the program for 
other ablative composites as their kinetic constants become 
available. 

The data shown in Tables 5-6 and 5-7 were used in this 
research to compute the degradation rate of the two compo- 
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sites. The decomposition rate computed from these data 
coupled to the char kinetic analysis (for the phenolic com- 
posite only) and the equilibrium analysis (for both compo- 
sites) constituted the novel approach in this research. The 
results of this approach are presented in the next chapter. 

Numerical Difficulties Encountered When Integrating The 
Chemical Kinetics Equations 

Numerical difficulties were encountered when the energy 
and species continuity equations for finite rate chemistry 
were simultaneously solved using a fourth order Runge-Kutta. 
The step size of the Runge-Kutta analysis was determined by 
the kinetics of the reacting flow. To preserve stability 
it was necessary , in many instances , to maintain a small 
step size, as small as 10”® feet, in regions where the flow 
was undergoing very rapid-chemical changes . This condition 
created by the very rapid chemical reactions caused the in- 
tegration to proceed very slowly. This numerical condition 
is known as stiffness (1, 10-15). It occurs in many physical 
systems which give rise to ordinary differential equations. 

It occurs when the simultaneous relaxation of different 
components vary at greatly different time rates. Mathemati- 
cally it is observed when the eigenvalues of the differen- 
tial equations are widely separated (13-15) . A parasitic 
eigenvalue is one that is associated with the extraneous 
solution of the difference equation used to approximate the 
differential equation. This extraneous solution grows with 
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time, and will eventually dominate the numerical solution 
if it is not suppressed. The larger the parasitic eigen- 
value, the more difficult is the suppression of the spurious 
solution. When the solution to the differential equation 
does not contain an exponential function the positive eigen- 
value is called parasitic saddle and the negative, parasitic 
stiff. The integration interval is, thus, determined by the 
fastest rate, and the region of integration is determined by 
the slowest rate (11) . 

Physical Concept of Stiffness : The physical interpre- 

tation of stiffness in a non-equilibrium reacting flow 
system is that the species are created and destroyed in a 
time scale that is orders of magnitude smaller than the time 
scale of the fluid particles. When the ratio of flow time 
to reaction time increases as the reaction is moving towards 
equilibrium, the parasitic eigenvalues increase, causing the 
equations to become stiff (21) . 

Mathematical Concept of Stiffness (20) : Consider A to 
be a nonsingular m x m matrix with eigenvalue i=l...m. 

R e ^ i > ^ ' . If max (R^A^) /min (R^A^) is a large 

number then the linear system of differential equations: 


x = Ax 


(5-11) 
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is said to be stiff. For nonlinear system a similar argu- 
ment applies. 

Most standard methods of numerical integration are not 
effective for stiff systems. Their deficiency lies in their 
inability to approximate the true value of the stiff compo- 
nent. When this is achieved (as for example in the Pade 
approximation used by Magnus and Schectner (34)), the numeri- 
cal difficulty is bypassed. For most standard methods, as 
is the case with the Runge-Kutta, the numerical integration 
is not effective for stiff systems. They require a stabi- 
lity condition of the form: 

max |h.AjJ < const (5-12) 


which makes h smaller as A^ increases. When the eigenvalues 
of A have negative real parts, the condition 5.8 prevails 
for the duration of a calculation, even though the solution 
has essentially no dependence on theA^ (for which I Rg A^| 
is large), except for a small initial interval. 

To circumvent the restriction 5.8 for stiff system, 
different integration schemes have recently been published 
in the literature (1, 11, 12, 13, 14,, 15, 17, 18, 19, 20, 

21, 22, 23, 24, 25, 26, 31, 32, 33, 34, 35). Numerical 
experiments with these schemes have shown them to be more 
effective, the stiff er the system (30) . We shall examine 
some of them in the following paragraphs. 


■» 
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Numerical Methods for Integrating Coupled Differential 
Equations With Varying Time Cons tan t¥ 

In problems of high temperature flows, the chemical 
rate equations sometimes involve both fast and slow reac- 
tions in some regions of integration. Integrating the che- 
mical kinetic equations has been a major source of difficul- 
ty because the very fast reactions determine the integration 
interval causing the step size to become very small, and 
thus, increasing considerably computing time. Three 
procedures to integrate these equations have been used to 
solve this type of problem. These are the explicit series 
solutions (11, 12, 13, 14, 16, 19, 21, 24, 27, 28, 29), the 
implicit series solution (1, 10, 13, 14, 15, 17, 18, 19, 21, 
22, 23, 25, 26, 30), and Moretti's (31) linearization scheme 
which is also known as the subdomain method, or the method 
of rational approximation (19) . One example of an explicit 
series solution is the Runge-Kutta analysis which has a 
small truncation error. As with all explicit techniques, 
this method has the advantage of being simple and straight- 
forward. It requires only first derivatives and knowledge 
of the value of the independent variable at the previous 
step. It is a marching procedure which gives the value of 
the independent variable as the integration proceeds. The 
implicit methods, on the other hand, are much more compli- 
cated to program since they require iterative solutions, 
and in most cases, they also require higher order deriva- 
tives, as in the method of Lomax (14) . However, they have 
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the major advantage of requiring a larger step size to 
maintain stability. Near equilibrium when the step size is 
reduced, the implicit methods lose their value in speed and 
computation. Recently Liu (21) investigated both implicit 
and explicit methods of solution for solving non-equilibrium 
boundary layer problems. He claimed that the Runge-Kutta 
method was as good as any to obtain converged solutions. 
However, computer time was slightly higher than the other 
implicit techniques he investigated. The limit of the step 
size in the implicit method is dependent upon the magnitude 
of error introduced by truncation of the series used to 
approximate the solution and not on stability consideration 
(19) . 

The third method of solution, rational approximation, 
was proposed by Moretti (31) . In it Moretti linearized 
the rate term in the species continuity equation. Moretti' s 
method presents the advantage of a smaller step size. How- 
ever, it is more complicated to program than the Runge-Kutta 
analysis . 

In conclusion, then, explicit techniques as the Runge- 
Kutta analysis are easier to program than the implicit 
techniques. However, they suffer from the disadvantage of 
requiring a smaller step size than implicit techniques. 

They also fail when near equilibrium conditions exist? the 
implicit techniques do permit small advances toward equili- 
brium "but at a price of excessive labor and computing time!' 
( 21 ). 
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Justification for Selecting the Fourth Order Runge-Kutta 
Analysis 

In the previous section several of the most important 
techniques for integrating stiff equations were mentioned. 

In this examination , it was indicated that implicit tech- 
niques were unconditionally stable, while explicit techniques 
were conditionally stable and required a small value of h to 
keep them within the stability region. Even with this draw- 
back, the fourth order Runge-Kutta was chosen as the integra- 
tion scheme for this research. For our problem the case of 
programming far outweighted the computational speed of a 
more difficult technique to implement numerically. Besides, 
the implicit technique for one dimensional flow would have 
been attractive if the integration had to be performed over 
a distance of several feet, but not over a distance of 1/4 
inch as in this research. Moreover, Lomax (14) in his 
review of integration schemes for integrating stiff equations 
stated that "if nothing special is known a priori about the 
differential equation, the standard fourth order Runge-Kutta 
is probably the best. It is self-starting, has low computing 
storage capabilities, is easy to program, has good accuracy 
0(h5), and... is more stable than any of the standard predic- 
tor corrector process". It is interesting to note also, that 
Liu (21) used a fourth order Runge-Kutta to solve the boun- 
dary layer equations after a detailed study of different 
integration schemes which were developed to handle stiff 
equations . 
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Selecting a Runge-Kutta Step Size to Maintain Stability 
Kliegel and Tyson (18) have shown a method for 
computing a step size h that would maintain the numerical 
integration for the Runge-Kutta technique to be stable. 
This technique consisted in defining a fluid mechanic cha- 
racteristic time, , and a chemical relaxation time, T c 
respectively as: 

T f - L (5-9) 

tt 


and 


T„ = Ci - Ci* (5-10) 

c ac L /dt 

where L is a length and u a velocity, is the actual 
concentration of species i, and is the concentration 
if the system were at equilibrium, and dC^/dt is the rate 
of reaction of said species. To maintain stability it was 
required that: 


Tf h < 5.6 (5-11) 

T c 

This is a simple and straightforward method for 
checking stability. However, it was estimated that to com 
pute at every step the equilibrium composition of all the 
species for our system would more than offset any gain of 
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computing an optimum step size. Hence, a more rapid method 
of solving for a stable step size was required even if it 
did not provide the optimum solution of the step size. 

Two empirical techniques were developed to compute a 
step size that would give a stable solution. One was to 
keep a double material balance, one on the gas species and 
one on the solid species. As long as these material ba- 
lances checked to four significant figures, the solution 
was observed to be stable. This technique was rather suc- 
cessful as long as all species were present in larger than 
trace amounts, or when reaction rates would not abruptly 
change. Unfortunately, this technique did not work well 
for very small mass fluxes (0.005). 

The other technique developed was as follows: 

Consider , the mass flux of species i and AN ir the change 
in mass flux of i due to a chemical reaction. If, 


AN i < 0.1 (5-12) 

N. 

l 

the solution was observed to be stable. 

This empirical criteria of Equation (5-12) works as 
long as the kinetics are controlling the step size. However 
it does not work at lower temperatures (1500°F) where the 
fluid mechanics step is the controlling element in the 
solution. This problem was circumvented by specifying a 
maximum step size that maintained stability when the reac- 
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tions were very slow. The computer running time was on the 
average , 15 minutes on an IBM 360 Model 65. 

Summary 

The development of the non-equilibrium flow analysis 
has been presented. The success of this analysis lies in 
part in selecting the important chemical reactions in the 
temperature range of interest, which in this research is 
500°F to 5500°F. To select these chemical reactions a 
three step procedure was outlined. First, the composition 
of the pyrolysis gases had to be established as accurately 
as possible. This was necessary to determine all possible 
chemical reaction combinations based on the original pyro- 
lysis products. Once these reactions were compiled, a 
thermodynamic equilibrium study was performed by computing 
the extent of each reaction by the free energy minimization 
technique. This second step helped to determine those reac- 
tions which were not thermodynamically feasible in the 
temperature range 500 to 5500°F. The final step was to 
analyze the kinetic data for each of the chemical reactions 
left. This was done by computing the conversion of an equi- 
molal mixture of reactants flowing isothermally through 
0.25 inch reactor. 

The second part of this chapter dealt with the numeri- 
cal implementation to describe the rate of finite reaction 
process. It was shown that the rate of reaction of a 
specie could be represented by a simple phenomenological 
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expression of the form of Equation (5-2) . The implementa- 
tion of this equation for computer use is simple, however, 
its solution was the cause of many numerical problems. It 
was seen that these numerical problems, associated with 
reacting flow, have been attributed by many as a mathemati- 
cal phenomena known as stiffness (1, 10-23, 31). Three 
numerical methods currently in use to integrate the coupled 
energy and chemical kinetic (continuity) equations were 
briefly described. These were the Runge-Kutta, implicit 
procedures, and the subdomain method, or method of rational 
approximation . 

In the next chapter the results of the solution of 
the equations of change for decomposition in depth will be 
presented. 
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CHAPTER VI 


ANALYSIS OF THE IN-DEPTH RESPONSE 
OF ABLATIVE COMPOSITES 

The ablation phenomena of heat shield materials during 
reentry is comprised of a large number of complex chemical, 
thermal and physical processes. Three analyses were deve- 
loped to predict heat shield performance, i.e., the total 
energy absorbed by an ablator. These analyses differ in the 
method of calculating the chemical generation term. When the 
chemical generation term is zero, the analysis is called 
frozen; when calculated from equilibrium thermodynamic con- 
siderations it is called equilibrium; when calculated using 
reaction rate information it is called non-equilibrium or 
kinetics analysis. These analyses have been described in 
previous chapters. 

The analyses developed in this research are to describe 
the energy absorbed in both the decomposition and char zone 
of char forming ablators. The first two analyses, frozen and 
equilibrium, were used to compute the lower and upper bounds 
of the energy absorbed by a nylon-phenolic resin composite and 
a silicon elastomer composite. The non -equilibrium flow ana- 
lysis was used to compute the energy absorbed in the decompo- 
sition zone and in the char zone of a nylon-phenolic resin 
composite only. This latter analysis was not 
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extended to the silicon elastomers for two reasons: One was 

the lack of good kinetic data on the silicon-carbon reactions 
and the other was the absence of reliable data on the compo- 
nents and composition resulting from the degradation of these 
polymers. Accurate information on the species resulting from 
the degradation of the polymers is essential to develop an 
accurate and realistic kinetics analysis. That is, to pos- 
tulate all the important chemical reactions requires an 
intimate knowledge of the species that form during polymer 
degradation. This information was not available for the 
development of a silicone elastomer kinetics model. 

Results of the energy absorbed as predicted by the 
three analyses are presented in this chapter. Of the three 
analyses, equilibrium and non-equilibrium analyses are very 
close to each other in the prediction of the total energy 
absorbed, while frozen predicts a total energy absorbed 
which is approximately a factor of 10 lower than the other 
two analyses. 

Energy Absorption in an Ablator 

We have discussed in previous chapters the foundation 
for the development of the three methods of analyses used 
in this research which are frozen, equilibrium and non- 
equilibrium. In this chapter, we are going to compare the 
results of the three analyses. We shall begin by explaining 
the principal energy absorbing mechanism in an ablator, both 
in the virgin plastic and in the char, and, to what extent 


254 


the three methods of analysis possess or lack these energy 
absorbing mechanisms. In addition, we will give an approxi- 
mate temperature range at which these energy absorbing me- 
chanisms become important contributing factors to the total 
energy absorbed in the ablator. 

Energy Absorbing Mechanisms in the Virgin Plastic 

In Chapter V we presented the data for the decomposi- 
tion of phenolic-nylon by Sykes and Nelson (3) . In Figure 
5-3 we showed the thermogram of nylon-phenolic and phenolic 
microballoon decomposition with the total energy absorbed 
for each of these ablator components. Using Sykes et. al . 

(3) data, we calculated, for example, that at a surface 
recession velocity of 0.02 ft/sec, the gas mass flux gene- 
rated by the decomposition of the virgin plastic for the 
non-equilibrium case is 0.4231 lb/ft^-sec. From this data 
we also calculated the total heat absorbed by the decomposi- 
tion process itself and found it to be approximately 91 BTU/ 
ft 2 -sec at the above mentioned surface recession velocity. 
The total energy absorbed in the virgin plastic, however, 
was approximately 99 BTU/ft^-sec. The difference of 8 BTU/ 
ft 2 -sec is the heat absorbed by the sensible enthalpy of the 
plastic, and of the sensible enthalpy of the gas generated 
by the plastic decomposition. 

In conclusion, we can say that the predominant energy 
absorbing mechanism in the plastic zone is that due to the 
heat absorbed by the decomposition process of the plastic 
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composites, and less than 10 percent to the sensible 
enthalpy gain of the plastic and the gas. In our analysis 
we did not segregate the decomposition process into a 
melting step and the actual decomposition or breakdown of the 
polymer. Rather, we lumped both together. However, by 
analyzing the data in Figure 5-3 we can see that, as expec- 
ted, most of the heat absorbed is due to the breakdown of 
the polymer and 10-12 percent only is due to melting of the 
polymer. 

Energy Absorbing Mechanisms in the Char 

The energy absorbing mechanisms in the char are 
basically three. One is due to sensible enthalpy gain, the 
other is due to chemical reactions, and the third one is 
due to sublimation. In the case of the virgin plastic the 
energy absorbing mechanisms were independent of whether the 
analyses used were frozen, equilibrium or non- equilibrium 
since the treatment of the virgin plastic is independent 
of the method of analyses in the char zone. However, in the 
char zone for example, chemical reactions are a function of 
which method of analysis we used. Therefore, identification 
of the temperature range at which each mechanism becomes 
important is much more difficult. 

First Energy Absorbing Mechanism ; Of the three energy 
absorbing mechanisms, two are related to chemical or phase 
changes. The energy absorbing mechanism in the frozen flow 
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case is that associated with the heat absorbed by the gas 
(transpiration cooling) in conjunction with the heat absorbed 
by the solid carbon matrix (the char) . Figure 6-1 shows a 
plot of the total energy absorbed as calculated by frozen, 
equilibrium and non-equilibrium case. As we can see from 
this plot, the total amount of energy absorbed by the frozen 
flow (transpiration cooling and sensible enthalpy gain of 
the char) is 1500 BTU/ft 2 -sec, while for non-equilibrium and 
equilibrium, the energy absorbed is 18300 and 19300 BTU/ft 2 - 
sec, respectively. We can see from this plot that transpi- 
ration cooling plays a very small part in the total amount 
of energy absorbed in the char. Examination of the data of 
Figure 6-1 shows that sensible enthalpy gain of the gas and 
the solid matrix is an order of magnitude~lower than that 
absorbed by chemical reactions. 

Second Energy Absorbing Mechanism : The second energy 

absorbing mechanism is that due to the heat absorbed by 
chemical reactions. This is best illustrated by comparing 
the difference in energy absorbed between equilibrium and 
frozen analyses. Figure 6-2 illustrates this difference. 

This large difference is due to the highly endothermic che- 
mical reactions that take place in the char zone. If this 
were not the case, the difference would not be that pronoun- 
ced. Equilibrium analysis, however, is an ideal model of the 
reactions, and comparison of frozen and equilibrium analyses 
throws little light as to when reactions, in practice, become 
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Figure 6-1. A Comparison of the Rate of Heat Absorbed 
Frozen, Non-equilibrium and Equilibrium 
Analyses at a Pressure of 0.1 Atm. and a 
Total Mass Flux of 2.1 lb/ft -sec., 
(v=0.06 ft/sec). 
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Figure 6-2. A Comparison of the Rate of Heat Absorbed 
Between Frozen and Equilibrium as a 
Function of Ablator Temperature at 0.1 
atmosphere and a Total Mass Flux of 2.10 
lb/ft -sec. (v=0.06 ft/sec). 
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important energy absorbing mechanisms. However, a comparison 
of frozen and non-equilibrium analyses or kinetic analysis 
should show the difference. For example. Figure 6-3 com- 
pares the total energy absorbed between frozen and non-equi- 
librium at a total mass flux of 2.10 lb/ft 2 -sec. Figure 6-4 
compares the two analyses also, but at a mass flux of 0.7 
lb/ft 2 -sec. At the higher mass flux of 2.10, the decomposi- 
tion/char zone temperature interface is approximately 2900° 

F, where the gases start reacting quickly as is evident by 
the sharp difference in slope. At the lower mass flux of 0.7 
where the decomposition/char zone temperature interface is 
lower (r* 2250°F) , the gases also start reacting quickly. 

The studies of April and Pike et. al. (2) also showed that 
reactions become important energy absorbing mechanisms at 
temperatures above 2000°F. To illustrate this point we show 
Figures 6-5, 6-6, 6-7 and 6-8. Figure 6-5 compares the 
temperature profiles for equilibrium, non-equilibrium and 
frozen flow at an assumed front surface temperature of 1500° 

F and an imposed mass flux of 0.05 lb/ft 2 -sec at the back 
surface of the char. This figure shows that both the frozen 
and non-equilibrium curves, for all practical purposes, are 
indistinguishable. This is because at the low temperatures 
of 1500°F the gases remain essentially frozen as the tempera- 
ture -is not high enough for chemical reactions to begin. 
Figure 6-6, on the other hand, shows that for an assumed 
front surface temperature of 2000°F the non-equilibrium and 
frozen flow analyses differ. Figures 6-7 and 6-8 which are 
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Figure 6-3. A Comparison of the Rate of Heat. Absorbed Between 
Frozen and Non-equilibrium at a Pressure of 0.1 ? 
Atmospheres and a Total Mass Flux of 2.10 lb/ft - 
sec. (v=0.06 ft/sec). 




HEAT ABSORBED (BTU/FT -SEC) 


261 



TEMPERATURE °F 

Figure 6-4. A Comparison of the Rate of Heat Absorbed 

for Frozen, Equilibrium and Non-equilibrium 
Analyses at a Pressure of 0.1 Atmospheres 
and a Total Mass Flux of 0.7 lb/ft -sec. 
(v=0.02 ft/sec). 
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Figure 6-5. Temperature Profile for the Frozen, 

Equilibrium, and Non-equilibrium Flow 
of Pyrolysis Gases Through the Char 
Zone of a Nylon-Phenolic Resin Ablator Q 
at a Front Surface Temperature of 1500 F 
as Reported by April (1) and Pike et. al. 
( 2 ). 
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Figure 6-6. 

of Pyrolysis Gases Through the Char 
Zone of a Nylon-Phenolic Resin Ablator 0 
at a Front Surface Temperature of 2000 F 
as Reported by April (l) and Pike et. al . 
( 2 ). 


Temperature Profile for the Frozen, 
Equilibrium, and Non-equilibrium Flow 
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Figure 6-7. Temperature Profile for the Frozen, 

Equilibrium, and Non-equilibrium Flow 
of Pyrolysis Gases Through the Char 
Zone of a Nylon-Phenolic Resin Ablator 0 
At a Front Surface Temperature of 2500 F 
as Reported by April (l) and Pike et. al. 
( 2 ). 
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Figure 6-8. Temperature Profile for the Frozen, 

Equilibrium, and Non-equilibrium Flow 
of Pyrolysis Gases Through the Char 
Zone of a Nylon-Phenolic Resin Ablator 0 
at a Front Surface Temperature of 3000 F 
as Reported by April (1) and Pike et. al. 
( 2 ). — — 
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plots for 2500 and 3000°F front surface temperature, respec- 
tively, further illustrate this point. Therefore, we can 

state that kinetic reactions become important energy ab- 

1 o 

sorbing mechanisms above 2000 F. 

Third Energy Absorbing Mechanisms ; In Figure 6-9 a 
comparison of the total energy absorbed versus temperature 
for equilibrium and non-equilibrium with the behavior of the 
carbon/ gas ratio superimposed, is shown. As shown in this 
figure, the carbon/gas ratio for both analyses remain 
practically constant up to a temperature of 4800°F, At this 
temperature the ratios begin to decrease, but with equili- 
brium decreasing at a faster rate. The reduction in carbon/ 
gas ratio is due in part to a mechanism known as sublima- 
tion, where the solid carbon becomes a gas. It is this 
mechanism (i.e., sublimation) which accounts for the di- 
vergence of the two analyses at the higher energy spectrum. 

In conclusion, we have shown that there are three 
principal energy absorbing mechanisms in the char: trans- 

piration cooling which is practically independent of 
temperature, except for the temperature effect on heat 

capacity; a reaction kinetic regime which begins at tem- 

o 

peratures above 2000 F; and a sublimation region which 

o 

begins at temperatures of about 4-800 F and above. 

A Comparison of the Decomposition/ Char Zone Boundary Condi- 
tions Between Equilibrium and Non-Equilibrium Analyses. ~ 


As we saw in Figure 6-1, the decomposition-char zone 
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Figure 6-9. A Comparison of Carbon/Gas Ratio for the 
Equilibrium and Non-equilibrium Analyses 
Illustrates the Effect of thg Divergence 
on the Two Curves Above 4800 F. 
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interface temperatures are different. For non-equilibrium 
this temperature is determined when the density of the de- 
grading plastic is equal to the density of the char. For 
phenolic-nylon the density of the virqin, undegraded polymer 
is approximately 35 lb/ft^. As the polymer is heated the 
virgin plastic begins to degrade and looses mass (as a gas) 
but not volume. Therefore, the density of the polymer be- 
gins to decrease as we integrate forward towards the char 
front. 

There is a point at which the polymer density 
becomes equal to the char density, and we call this the 
decomposition/char zone interface. For phenolic-nylon this 
char density has been experimentally determined to be 13 lb/ 
ft^. Therefore when the polymer density reaches 13 lb/ft^, 
the program enters the char zone and activates the kinetic 
model . 

In the case of the equilibrium analysis we cannot 
use the density of 13 lb/ft^ as the criteria for defining 
the decomposition/char zone interface. If we were to use 
this criteria we would find a discontinuity at the boundary 
between the two reqions. This discontinuity would be in the 
mass flux of the gas predicted by the polymer decomposition 
kinetics and that predicted by the equilibrium analysis. 

To avoid this discontinuity we define the interface at that 
point where the mass flux from the decomposition of the 
polymer equals that required by chemical equilibrium. Since 
the equilibrium composition of a mixture is only a function 
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of temperature, pressure and elemental composition, we can 
determine a priori, for phenolic-nylon, what the equilibrium 
mass flux has to be as a function of temperature for every 
surface recession velocity. Having done this, the computer 
program checks at every integration step whether or not the 
mass flux generated by polymer degradation is equal to that 
required by equilibrium considerations. When these condi- 
tions are met for equilibrium analysis, we define this as 
the decomposition/char zone interface. When the integration 
procedure enters the char, the gases are assumed to be in 
chemical equilibrium and the program continues to integrate 
the equations of change under this assumption. 

It is interesting to note at this point that be- 
cause equilibrium tends to overpredict the carbon/gas ratio 
as compared to actual experiments, the mass flux of gas 
formed by the degradation of the polymer required to match 
equilibrium conditions, is small. Therefore, as shown in 
Table 6-1 (P=0.1 atm.) we see that at the end of the decom- 
position zone or beginning of the char zone, the gas mass 
flux for equilibrium is smaller than for non-equilibrium. 

As we see the gas mass flux for non-equilibrium is approxima- 
tely 61.6 percent greater than for equilibrium for a total 
mass flux pv, of 0.7 lb/ft 2 -sec., (or a surface recession 
velocitv, v, of 0.02 ft/sec) . At a total mass flux of 2.10 
(v=0.06 ft/sec) the gas mass flux for non-equilibrium is 
61.4 percent qreater. Table 6-2 is similar to 6-1 but the 
conditions are those for a pressure of 1 atmosnhere. At a 



TABLE 6-1 : Comparison of Equilibrium and Non-Equilibrium Analysis 

Gas Mass Flux at the Back Surface of the Char, at Various 
Surface Recession Velocities* and at a Pressure of 0.1 
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surface recession velocity of 0.02 ft/sec non-equilibrium 
gas mass flux is 63.2 percent greater than for equilibrium. 
Therefore, the gas mass flux at the back surface of the char 
is always greater for non-equilibrium than it is for equili- 
brium. It is important that we keep in mind this fact, 
because it bears heavily on explaining the behavior of the 
equilibrium and non-equilibrium heat curves plots shown in 
this chapter. 

Table 6-3 and 6-4 show the decomposition/char zone 
temperature interface at various surface recession veloci- 
ties for 0.1 and 1.0 atmospheres respectively. We see in 
both of these tables that the temperature of the interface 
is lower for equilibrium than for non-equilibrium. This is 
expected since as we explained before, the extent of polymer 
degradation required to meet the gas mass flux at the equi- 
librium boundary conditions is less than for non-equilibrium. 
For example, at a surface recession velocity of 0.02 ft/sec. 
the temperature required to meet the equilibrium boundary 
conditions is 1386. 9°F as is shown in Table 6-3. This is 
929°F less than for non-equilibrium. In comparing the two 
tables we can also notice that the temperature at the inter- 
face is not affected by pressure. This is expected because 
the depolymerization kinetics are not a function of pressure, 
and as we said before, the plastic will keep degrading until 
it reaches a density of 13 lb/ft , which is the density of 
the char. On the other hand, pressure slightly affects the 
interface temperature in all the cases examined for equili- 



273 


n 

•H X 

o) m 
tx 

f-« •» 

m <d 
c c 
< o 

N 

x 

ia 


E 

3 

•X 


X 

o 

<1) 

X 

3 

0) 

0) 

<D 

X 

cu 


« 

C 

X U 

0 

•X *0 

2 

rl 0) C 


•X X3 03 


3 X 


f 01 


W X <D 


1 0 -X 


c: x 

01 

0 QJ -X 

•H 

z: o o 

01 

(0 0 

>i 


i— 1 

C5 Wt <U 

03 

Q3> 


U) 


e e 


3«o 

£ 

•X 0 -X 

3 

X 03 01 

•X 

X ffl 01 

X 

•X 0) 

X 

H 51 U- 

•X 

•X X 0) 

fX 

3 x a 

•X 

3* 

3 

W X <l) • 

cr 

osuai 

w 

X 03 X 


o ai x a> 


01 x X 


c x 3 a, 


0 3 w 01 


01 X O 

01 

-X 03 01 E 

0) 

X X 3 x 

•X 

03 0) O 03 

X 

Qj Qj “X 

•X 

E E X rX 

o 

0 0) 03 • 

•X 

U H > o 

f — I 


0) 


> 

• • 

c 

cn 

o 

l 

•X 

vo 

01 


01 


<D 

M 

o 

X 

0) 

a 

m 

< 



01 


o 


03 


X 


X 


<3 


01 

•H 

0) 

>1 

r— 1 

03 


•H 

X 

XI 
•X Cli 
rX 3 
•H 
3 

cr 
w 

I 


0 

QJ 

01 


o\ 

VO 

in 

rX 

r- 

• 

• 

• 

• 

• 

in 

HJ* 


ov 

m* 

iX 

r-l 

VO 

00 

cn 

<n 

in 

VO 

I" 

00 

CM 

(N 

CM 

CM 

CM 


OV 

r-~ 

ov 

CM 

a\ 

• 

• 

• 

• 

• 

VO 

o 

Cl 

r» 

o 

00 

m 

OV 

m 

r- 

m 

M* 

M* 

in 

in 

H 

«X 

rX 

rX 

i— < 


CM 

cn 


m 

VO 

O 

o 

o 

o 

o 

d 

d 

d 

d 

d 



274 


3 <4-1 

•r< 4J O 
CO (0 


>« 3 

rH ^ 4-1 

<0 <U 3 
3 3 to 

< O CO 

N <D 

g n 

3 Jh 04 
■H 3 
H .3 3 
XI CJ 

•r< T} 

H J) C 
■H £ ID 
3 -3 
tr to 
K*w ® 

1 O *H 
3 +» 

O 3-H 
ZOO 
3 O 
tJ<WH 
3 3 0 
3 3 > 

CO 

g 3 
3 O 
■rl O -H 
3 3 3 
XJ ffi m 

•rH a) 

H dl u 

•w si a) 

3 -*j os 
o - 

W +> 3 • 
3 0 3 
4-) (T 3 
0 3 4-13 
3 3 X! 
CH3ft 
O 3 CO 3 
3 -P o 
t4 3 3 E 
3 3 3 -P 
3 3 O 3 
& Q -H 
g £ 3 O 
0 3 3* 
O Eh > rH 


I 

\D 

w 

31 

w 

< 

Eh 





275 


brium; the difference varying from a high of 2.7°F at 
v=0.02 ft/sec to a low of 0.09°F at v=0.06 ft/sec (see Tables 
6-3 and 6-4) . It should also be noted that the interface 
temperature monotonically increases with surface recession 
velocity for both equilibrium and non-equilibrium. This 
should be expected since the higher the surface recession 
velocity (or the total mass flux), that is, the more severe 
the conditions that we are trying to simulate are, then the 
steeper the temperature profile should be inside the ablator. 

Table 6-5 and 6-6 show the chemical composition of 
the gases at the char back end for 0.1 and 1 atmospheres 
respectivelv. Each shows in addition one species composi- 
tion at two surface recession velocities. As we can see in 
Table 6-5, the species composition for equilibrium varies 
depending on whether the surface recession velocity v, is 
0.02 or 0.06 ft/sec. This is not because the compositions 
are a direct function of surface recession velocity, but 
because the compositions are a function of back surface tem- 
perature, which in turn is a function of surface recession 
velocity. (See Tables 6-3 and 6-4 which show how surface 
recession velocities affect decomposition/char zone tempera- 
ture interface) . As we examine Table 6-5 and 6-6, we can 
notice that the equilibrium compositions are slightly 
different at the two surface recession velocities, but yet 
do not vary for non-equilibrium. The same is true when we 
compare them at both 0.1 and 1.0 atmospheres respectively. 
Unfortunately, this should not be so, but as we have said 
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4 

TABLE 6-5: A Comparison of Species Composition at the Back 

Surface of the Char for Equilibrium and Non-Equilibrium 
Analyses at a Surface Recession Velocity of 0.02 and 
0.06 ft/sec and a Pressure of 0.1 Atmosphere. 


Equilibrium Non-Equilibrium 



(mole 

fraction) 

(mole 

fraction) 


FT/SEC: 
v= 0 . 02 

v=0 . 06 

FT/SEC: 
v= 0 . 02 

( 

v=0 .06 

CH, 





CH 3 

- 

- 

- 

- 

CH. 

0.0041 

0.0008 

0.0804 

0.0804 

C 2 H 

- 

- 

- 

- 

C 3 H 

- 

— 

- 

- 

c 4 h 

— 

— 

— 

— 

C 2»2 

— 

— 

0.0480 

0.0480 

C 2 H 4 

- 

- 

0.0450 

0.0450 

C 2 H 6 

- 

- 

0.0070 

0.0070 

C 6 »6 

- 

— 

0.0107 

0.0107 

c 6 h 6 o 

— 

— 

0.0794 

0.0794 

CN 

- 

- 

- 

- 

CO 

0.2105 

0.2153 

0.0480 

0.0480 

CO 2 

0.0015 

0.0001 

0.0338 

0.0338 

H 

— 

— 

— 

— 

«2 

0.7379 

0.7419 

0.4891 

0.4891 

HoO 

0.0041 

0.0003 

0.1014 

0.1014 

OH 

- 

- 

- 

- 

M 2 

0.0419 

0.0416 

0.0572 

0.0572 

mh 3 

- 

- 

- 

- 

HCN 

- 

- 

- 

- 

0 

- 

- 

- 

- 

C 3 

- 

- 

— 

- 

C (S) * 

1.227 

1.215 

1.062 

1.062 


*The symbol C( s ) represents the moles of carbon per mole of 
gas . 
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TABLE 6-6: A Comparison of Species Composition at the 

Back Surface of the Char for Equilibrium and Non- 
Equilibrium Analyses at a Surface Recession 
Velocity of 0.02 and 0.06 ft/sec and a Pressure 
of 1.0 Atmosphere. 



Equilibrium 
(mole fraction) 

Non-Equilibrium 
(mole fraction) 


FT/SEC : 
v= 0 . 02 

v=0 . 06 

FT/SEC : 
v= 0 . 02 

v=0 . 06 

ch 7 

_ 



«• 

ch 3 

- 

- 

- 

- 

ch 4 

0.0366 

0.0151 

0.0804 

0.0804 

c 2 h 

- 

- 

- 

- 

c 3 h 

— 

— 

— 

- 

C 4 H 

- 

- 

— 

- 

C 2 H 2 

- 

- 

0.0480 

0.0480 

c 2 h 4 

- 

— 

0.0450 

0.0450 

C 2 h 6 

- 

- 

0.0070 

0.0070 

C 6 H 6 

- 

— 

0.0107 

0.0107 

c 6 h 6 0 

- 

- 

0.0794 

0.0794 

CN 

- 

- 

- 

— 

CO 

0.1775 

0.2071 

0.0480 

0.0480 

co 2 

0.0110 

0.0025 

0.0338 

0.0338 

H 

— 

— 

— 

— 

H 2 

0.6964 

0.7234 

0.4891 

0.4891 

h 2 o 

0.0335 

0.0094 

0.1014 

0.1014 

OH 

- 

- 

- 

- 

n 2 

0.0449 

0.0425 

0.0572 

0.0572 

nh 3 

0.0001 

- 

- 

- 

HCN 

- 

- 

- 

- 

0 

— 

— 

— 

— 

c 3 

- 

— 

- 

- 

c (s)* 

1.320 

1.24 

1.062 

1.062 

* The 
gas. 

symbol C( s j 

represents 

the moles of carbon 

per mole 





278 


before, there is no known method available to predict non- 
equilibrium composition of degrading gases directly from the 
polymer degradation. Therefore, as a compromise we had to 
use available experimental data and come up with the best 
estimate of the species composition at the back surface of 
the char zone. Therefore, for all non-equilibrium cases run, 
varying the surface recession velocity and the pressure, we 
used the same chemical composition. 

A Comparison of the Total Energy Absorbed for Equilibrium , 
Non-Equilibrium and Frozen Flow Analysis for Phenolic-Nylon 

In developing the three methods of analysis, frozen, 
equilibrium and non-equilibrium, the purpose of this research 
was to determine the effects that each method of analysis 
had on the predicted value of the total energy absorbed in 
an ablator. Before we get into the details of analyzing 
the cases examined in this research, it is worthwhile to 
briefly summarize how the approach of our analysis evolved 
and how it differs from the initial research of April (l)and 
Pike et. al. (2) . The initial research that was done invol- 
ved looking at the char zone only. In this initial stage of 
development a mass flux was assumed at the back surface of 
the char, and in addition a front surface temperature of the 
char was assumed as a parameter (see Figures 6-5 through 
6-8) . The back surface of the char was always assumed to be 
at 500°F. Therefore, in comparing frozen, equilibrium and 
non-equilibrium calculations, the boundary conditions at the 
back end of the char were always the same in terms of 
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temperature and mass flow. 

In this research we went a step further and coupled the 
plastic zone to the decomposition zone. The result of this 
coupling makes the analysis of the heat curves more complex. 
Basically , because this coupling, as we have seen in Tables 
6-1 through 6-6, has resulted in boundary conditions which 
are not the same for equilibrium and non-equilibrium both 
in terms of temperature and mass flux. As a result when we 
look at the total energy absorbed we have to take these 
differences into consideration in analyzing the values of the 
total energy absorbed. 

In Figures 6-10 through 6-16 we present the comparison 
of the total energy absorbed for the three methods of ana- 
lysis at two pressure levels. As we examine these curves we 
notice that the non-equilibrium curves cross the equilibrium 
curves for all cases examined, except one. The reason for 
this crossover is obvious when we look at the gas mass flux 
at the back surface of the char. As we have already shown, 
the gas mass flux for non-equilibrium is about 60 percent 
greater than for equilibrium. The crossover occurs because 
at temperatures above 2000°F the total rate of heat absorbed 
by the non-equilibrium gases is initially greater than for 
equilibrium. This is not because the non-equilibrium analy- 
sis chemical reactions are more endothermic than those for 
chemical equilibrium, but because the higher gas mass flux 
of the non-equilibrium analysis coupled with the heat ab- 
sorbed by the kinetic reactions is sufficient to overcome the 
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500 1000 2000 3000 4000 5000 6000 

TEMPERATURE °F 


Figure 6-10. A Comparison of the Rate of Heat Absorbed 

for Frozen, Equilibrium and Non-equilibrium 
Analyses at a Pressure of 0.1 atm and a Mass 
Flux of 0.7 lb/ft z -sec (v=0.02 ft/sec). 
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TEMPERATURE °F 

Figure 6-11. A Comparison of the Rate of Heat Absorbed 

for Frozen. Equilibrium and Non-equilibrium 
Analyses at a Pressure of 0.1 atm and a 
Total Mass Flux of 1.05 lb/ft^/sec). 
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Figure 6-12. A Comparison of the Rate of Heat Absorbed 

for Frozen, Equilibrium and Non-equilibrium 
Analyses at a Pressure of O.l^atm and a 
Total Mass Flux of 1.40 lb/ft -sec, 

(v=0.04 ft/sec). 
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IOOO 2000 3000 4000 5000 6000 

TEMPERATURE °F 

Figure 6-13. A Comparison of the Rate of Heat Absorbed 

for Frozen, Equilibrium and Non-equilibrium 
Analyses at a Pressure of 1.0 atm and a 
Total Mass Flux of 0.7 lb/ft -sec, 

(v=0.02 ft/sec). 
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Figure 6-14. A Comparison of the Rate of Heat Absorbed 

for Frozen, Equilibrium and Non-Equilibrium 
Analyses at a Pressure of 1.0 atm and a 
Total Mass Flux of 1.05 lb/ft^-sec, 

(v=0.03 ft/sec). 
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Figure 6-15. A Comparison of the Rate of Heat Absorbed 

for Frozen, Equilibrium and Non-equilibrium 
Analyses at a Pressure of 1,0 atm and a 
Total Mass Flux of 1.4 lb/ft 2 -sec, 

(v=0.04 ft/sec). 
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Figure 6-16. A Comparison of the Rate of Heat Absorbed 

for Frozen, Equilibrium and Non-equilibrium ' 
Analyses at a Pressure of 1.0 atm and a 
Total Mass Flux of 2.1 lb/ft^-sec, 

(v=0.06 ft/sec). 




287 


much lower mass flux rate of chemical equilibrium. As the 
temperature increases, however, the carbon in the chemical 
equilibrium model begins to react into the gas phase, in- 
creasing the equilibrium gas mass flux. As this same mass 
flux increases, the equilibrium heat absorbed finally over- 
takes the non-equilibrium curve and surpasses it. 

A Comparison of the Total Heat Absorbed for Frozen, Equili- 
brium and Non-Equilibrium at 0.1 Atmosphere 

The results for frozen, equilibrium and non-equilibrium 
analyses have been calculated at various surface recession 
velocities and at two pressure levels. In this section we 
will discuss the results at 0.1 atmosphere. 

Figure 6-10 shows a plot of the total energy absorbed 

as calculated by the three methods of analysis at a surface 

recession velocity, v, of 0.02 ft/sec. This velocity 

corresponds to a total mass flux of 0.7 lb/ft -sec. It 

should be noted again that the total mass flux is the product 

3 

of the density of the virgin material (35 lb/ft for phenolic- 
nylon) and the surface recession velocity. The arrows shown 
in this figure, and subsequent figures, denote the end of the 
decomposition zone or beginning of the char zone. The reader 
is referred to Table 6-3 for the exact temperatures of decom- 
position/char zone interface. Figures 6-11, 6-12 and 6-3 
(which has been shown earlier) correspond to total mass flu- 
xes, Pv, of 1.05(v=0.03 ft/sec), 1.40(v=0.04 ft/sec) and 
2.10(v=0.06 ft/sec) lb/ft^-sec. We should notice that all 
the equilibrium curves shown in these figures, end at a tem- 
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perature of 5500°F. This temperature was not selected ar- 
bitrarily, rather it is the temperature at which all the 
solid carbon, either by reacting with other gases or by 
sublimation, disappear. The disappearance of all the solid 
carbon obviously denotes the physical end of the char zone. 
Table 6-7 shows a summary of the equilibrium flow results. 
Table 6-8, however, compares non-equilibrium and equilibrium 
flow results for the various surface recession velocities 
already mentioned. The reader will notice that the case 
for v=0.01 ft/sec was not run for the non-equilibrium case. 

The reason for this was the excessive amount of computer 
time required to complete this case. We will expand on this 
at the end of the chapter. 

2 

For the cases where ^v=0.7, 1.05 and 1.40 lb/ft -sec, 
there is a crossover between the equilibrium and non-equili- 
brium curves. We have already explained the reason for this 
crossover. This crossover however, does not occur for the 
case where total mass flux is 2.10 lb/ft^-sec (v=0.06 ft/sec). 
If we go back to Table 6-3 we will notice that as the surface 
recession velocity increases the difference between the equi- 
librium and non-equilibrium decomposition/char zone tempera- 
ture interface also increases; being smallest at v=0.02 ft/ 
sec, with a difference of 929°F, and largest at v=0.06 ft/sec, 
with a difference of 1324°F. This very large difference in 
temperature allows the equilibrium heat curve a head start. 

The non-equilibrium curve can never surpass it because by the 
time it approaches the equilibrium curve at about 4200°F, the 



TABLE 6-7: Comparison of the Total Energy Absorbed at Various Surface 

Recession Velocities: Equilibrium Analysis (P=0.1 Atm; 

Nylon-Phenolic Composites) . 
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carbon in the chemical equilibrium model begins to react 
into the gas phase increasing both the gas flux and the rate 
of heat absorbed. 

A Comparison of the Total Heat Absorbed for Frozen, Equili- 
brium and Non-Equilibrium Analyses at One Atmosphere 

In the previous section we discussed the results of 
frozen, equilibrium and non-equilibrium analyses for a pres- 
sure of 0.1 atmosphere. In this section we will discuss 
similar results but for a pressure of 1 atmosphere. 

Figure 6-13 shows a plot of the total heat absorbed by 
frozen, equilibrium and non-equilibrium for a surface reces- 
sion velocity of 0.02 ft/sec (fv=0 . 7 lb/ft 2 -sec) . The total 

amount of heat absorbed predicted by the equilibrium analy- 

2 

sis is 5599 BTU/ft -sec, that for non-equilibrium is 4800 
BTU/ft 2 -sec and for frozen, 625 BTU/ft 2 -sec. It should be 
noted that for 1.0 atmosphere the equilibrium analysis 
predicts a front surface temperature of 5800°F. This is the 
temperature at which the solid carbon concentration approaches 
zero. This temperature is 300°F higher than the case of 0.1 
atmosphere. One would expect this because when the pressure 
increases by a factor of ten, the energy required for the 
carbon to enter the gas phase has to be greater. Putting it 
another way, the higher the pressure, the higher the front 
surface temperature required to get the carbon into the gas 
phase. This is because the higher pressure tends to keep the 
carbon in the solid phase. Figures 6-14, 6-15 and 6-16 show 
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plots as those shown in Figure 6-13, but these are for 
surface recession velocities of 0.03, 0.04 and 0.06 ft/sec 
respectively. We have summarized the results of the total 
energy absorbed by equilibrium in Table 6-9. In Table 6-10 
we have summarized the results of the non-equilibrium ana- 
lysis and compared it to those of equilibrium. 

Effect of Surface Recession Velocity on Heat Absorbed on the 
Phenolic-Nylon Resin; Parameter Study 

In Table 6-8, the tabulated results of the total energy 
absorbed at a number of surface recession velocities are 
shown for both the equilibrium and non-equilibrium analyses. 
The pressure at which these results are tabulated is 0.1 
atmosphere. In Table 6-10 similar results are tabulated but 
for a pressure of 1 atmosphere. 

First we shall analyze the results at a pressure of 0.1 
for equilibrium. Figure 6-17 shows these results graphically. 
As expected, the higher the surface recession velocity, the 
larger the amount of heat absorbed. It should also be noted 
that the curves get closer to each other as the surface re- 
cession velocity increases. 

If we take the energy absorbed at v=0.01 ft/sec as the 
equilibrium base case shown in Table 6-6, and compare it to 
the cases for v=0.02, 0.03, 0.04, 0.05 and 0.06 ft/sec, 
there is almost a one to one correspondence between the ratio 
of the velocities and that of the total energy absorbed. So 
for 2, 3, 4, 5 and 6 times the surface recession velocity of 
the base case, we find that the total energy absorbed is 1.99, 




TABLE 6-9: Comparison of the Total Energy Absorbed at Various Surface 

Recession Velocities: Equilibrium Analysis (P=l. atm; 

Nylon-Phenolic Composites) . 
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Figure 6-17. A Comparison of the Total Energy Absorbed at Various 

Surface Recession Velocities: Equilibrium Analysis 

(P=0.1 atm; Nylon-Phenolic Composites). 
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2.98, 3.97, 4.95 and 5.94 times that of the base case 
respectively. One would expect it to be exactly 2, 3, 4, 
etc., times the base, that is, if we were analyzing the char 
zone only. After all, equilibrium reactions are not a func- 
tion of surface velocity, so doubling the mass through the 
char would exactly double the total amount of energy 
absorbed. However, the reason the ratios are not exactly 
the same can be found by looking at the decomposition/char 
zone interface temperature at each of these surface recession 
velocities (see Table 6-3) . We see for example that at 
v=0.02 ft/sec the interface temperature is 1387°F, at v=0.03 
ft/sec it is 1451°F, at v=0.04 ft/sec it is 1500°F, and so 
on. Therefore, the lower the surface recession velocity, 
the lower the back surface temperature of the char and the 
sooner the equilibrium gases will begin to react in the char. 
We can see, for example, that at vs0.02 ft/sec the gases are 
in equilibrium at 1387°F and begin to absorb heat at that 
temperature. At v=0.03 ft/sec, on the other hand, the gases 
will not reach equilibrium conditions until they reach a 
temperature of 1451°F. 

In Table 6-8, we show also the tabulated results of the 
total energy absorbed for non-equilibrium. In addition, the 
results are plotted in Figure 6-18. 

The non-equilibrium case shows more vividly the effect 
that the surface recession velocities have on the temperature 
boundary conditions and therefore on the shape of the heat 
curve. It also shows that an increase in the surface reces- 
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Figure 6-18. A Comparison of the Total Energy Absorbed at 
Various Surface Recession Velocities. 
Non-equilibrium Analysis (P=0.1 atm, Phe- 
nolic-Nylon Composites). 
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sion velocity results also in an increase in the temperature 
of the decomposition/char zone interface. This causes the 
heat curves to cross over one another close to the interface. 
For example, from Figure 6-18 we can see that the heat ab- 
sorbed at a surface recession velocity of 0.02 ft/sec is 
greater than that of the other surface recession velocities 
plotted between the temperatures of 2500°F to 2600°F. We 
can also see, for example, that at a surface recession velo- 
city of 0.03 ft/sec the accumulated heat absorbed is greater 
than that of 0.06 ft/sec, between the temperatures of 2600° 

F and 2900°F. The reason for this crossover is that at the 
lower mass velocities the gases begin to react earlier than 
at the higher mass flux. 

If we compare the energy absorbed at v=0.02 ft/sec with 
that at 0.03 ft/sec we see that the total heat absorbed at 
0.03 ft/sec is about 1.5 times that of 0.02 ft/sec. Since we 
do not have a case for v=0.01 ft/sec, for the reason that we 
have already mentioned, we will use v=0.02 ft/sec as the 
base case for comparison. As with equilibrium there is 
almost a one to one correspondence between an increase in 
surface recession velocity and total heat absorbed. We can 
see, therefore, that at 1.5, 2, and 3 times the surface re- 
cession velocity of 0.02 ft/sec, the total quantity of ener- 
gy absorbed is 1.49, 1.97 and 2.84 for 0.03, 0.04 and 0.06 
ft/sec respectively. In the equilibrium case we explained 
that the reason the correspondence of the ratios of energy 
absorbed were not exact integers of the surface recession 
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velocity (or for that matter the mass flux , fv) was 
because of the differences encountered in the boundary condi- 
tions . For non-equilibrium we can also apply the same 
reasoning. However, it should be noted that for non-equili- 
brium it is not necessary to get a close one-to-one corres- 
pondence. Basically because one can argue that the higher 
the surface recession velocity, i.e., the higher the mass 
flux through the char, the lower the residence time in the 
char and hence, the lower the time for reaction. Apparently, 
the reduction in residence time experienced at 0.06 ft/sec, 
for example, has had a small effect on the extent of the 
kinetic reactions. 

Effect of Pressure on Heat Absorbed by the Phenolic-Nylon 
Resin 

A parameter study was conducted on pressure to analyze 
the effect on the total heat absorbed. The reason for this 
study is that the reentry pressure is not constant and varies 
with the trajectory of the vehicle. Two pressure levels 
were selected to bracket this effect: these pressures were 

0.1 and 1.0 atmospheres respectively. The lower pressure 
level is the one encountered by the reentry vehicle at the 
higher altitudes while the 1.0 atmosphere level is the 
theoretical maximum and achieved only at the lower altitude 
in the trajectory. 

A plot of pressure effect for a surface recession velo- 
city of 0.06 ft/sec is shown in Figures 6-19 and 6-20. These 
plots are for the equilibrium and non-equilibrium analyses 




Figure 6-19. Effect of Pressure on the Total Heat Absorbed 

Equilibrium Analysis (v=0.06 ft/sec). 
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Figure 6-20 


Effect of Pressure on the Total Heat 
Absorbed.: Non-equilibrium Analysis 
(v=0.06 ft/sec). 
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respectively. 

For the equilibrium case shown in Figure 6-19, the 
heat curve for 1.0 atmosphere remains above that for 0.1 
atmosphere until it reaches a temperature of approximately 
3600°F, and then it falls below the heat curve for 0.1 
atmosphere. For the non-equilibrium case shown in Figure 
6-20 a similar behavior is manifested by the 1.0 atmosphere 
heat curve; although, the 1.0 atmosphere curve falls below 
the 0.1 atmosphere at a temperature of about 4600°F. Before 
the crossover for both the equilibrium and non-equilibrium 
case, the effect of pressure is small. 

Looking at Table 5-2, where the kinetic reactions are 
listed, one would expect that at P=0.1 atmosphere the heat 
curve would be higher since more of the kinetic reactions 
are favored by lower pressure than by higher pressure. 

There are some that are not affected at all. However, there 
are two reactions which are favored by higher pressure and 
these are the carbon-water and the carbon-carbon dioxide 
reactions. To be able to assess whether these reactions 
play an important part at temperatures lower than 4600°F we 
have to compare a plot of the species composition concentra- 
tion at both pressures. In Figures 6-21 and 6-22 a concentra- 
tion versus temperature plot is shown. In comparing the two 
plots we can see that the water concentration decreases more 
rapidly at 1.0 than at 0.1 atmosphere. In addition, it is 
also evident from the plot that the carbon-carbon dioxide 
reaction is taking place as is evident by the more rapid 
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TEMPERATURE °F 

Figure 6-21. Chemical Composition of the Phenolic-Nylon 
Pyrolysis Gases as Predicted by the Non- 
equilibrium Analysis at a Pressure of 1.0 
atm. and a Surface Recession Velocity of 
0.06 ft/sec. 
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Figure 6-22. Chemical Composition of the Phenolic-Nylon 

Pyrolysis Gases as Predicted by the Non- 
equilibrium Analysis at a Pressure of 0.1 
atm and a Surface Recession Velocity of 
0.06 ft/sec. 
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increases in carbon monoxide concentration. We see a similar 
behavior for the water and carbon dioxide concentration in 
the equilibrium calculations as is evident from the plots of 
Figure 6-23 and 6-24. Figure 6-23 and 6-24 are plots of 
equilibrium concentration versus temperature at pressures of 
1.0 and 0.1 atmospheres respectively. We want to caution 
the reader that the extrapolation of the non-equilibrium ar- 
gument to the equilibrium domain might be erroneous. In the 
calculations of equilibrium composition we used the free 
energy minimization technique which does not require the 
postulation of a reaction mechanism, and, as in any chemical 
equilibrium process, it is independent of the path. There- 
fore, to associate any particular mechanism to such complica- 
ted process may be dangerous. For the equilibrium case it 
should be sufficient to say that because of the free energy 
of the complex mixture, the endothermicity of the gas at 
a 1.0 atmosphere is greater below a temperature of 3600°F. 

What Figures 6-19 and 6-20 illustrate is that unless 
the front surface temperature is above 4000°F, the effect of 
pressure on total heat absorbed is essentially the same for 
0.1 and 1.0 atmosphere. 

Numerical Difficulties 

As we have mentioned before, the numerical integration 
technique used to solve the equations of change for frozen, 
equilibrium and non-equilibrium analyses was a fourth order 
Runge-Kutta. With the equilibrium and frozen analyses we 
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Figure 6-23. Equilibrium Composition of the Phenolic-Nylon 
Pyrolysis Gases in the Char at a Pressure of 
1.0 Atmospheres. 
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Figure 6-24. Equilibrium Composition of the Phenolic-Nylon 
Pyrolysis Gases in the Char, at a Pressure of 
0.1 Atmospheres. 
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did not experience any particular difficulty. However, with 
the non-equilibrium analysis we did. Basically, the reason 
for our difficulties were the stiff system of equations en- 
countered in the non-equilibrium analysis. Because of this 
stiff system the size of the integration step had to be very 
small to avoid a totally erroneous solution. As a result we 
used excessive amounts of computer time and had to limit the 
number of cases analyzed. In Figure 6-25 we show a plot of 
mass flux versus IBM/360 CPU time. As expected, at the 
higher mass flux the CPU time decreases in an exponential 
fashion. The reasons are basically two. One is that at the 
higher flow rates, the residence time is lower in the char. 
Secondly, as we showed earlier in the chapter, at the higher 
flow rates (mass flux) the temperature of the decomposition 
zone is higher and the thickness of the char is smaller. 
Therefore, the number of integration steps required to solve 
the equations of change are fewer, thus requiring a lesser 
amount of CPU time. 

Effect of Chemical Reaction Rate Data on the Non-Equilibrium 
Flow Calculations 

In order to study the sensitivity of the analysis to 
the chemistry model, a study was performed by taking each of 
the four reactions shown in Table 5-3 one at a time and using 
the reported sets of kinetic data discussed in Chapter V. 

The predicted values of the energy absorbed for each set of 
kinetic data was determined, and the conditions selected for 
the study were for a pressure of one atmosphere and a surface 
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Figure 6-25. Runge-Kutta Solution Time of the 
Non-equilibrium Analysis versus 
Mass Flux. 
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recession velocity of 1 x 10 ft/sec. For all the reactions 
in Table 5-3, except for the acetylene reaction, the varia- 
tions in the predicted values of the energy absorbed were 
less than one half of one percent. For the acetylene reac- 
tion the variation in energy absorbed was about 4 percent as 
discussed in Chapter V. These variations were not believed 
to be significant when one considers the wide discrepancy 
that exists in the kinetic data that is available in the 
literature. This close agreement is not fortuitous. But it 
is the product of a very critical study of the kinetic 
literature, of a very meticulous approach in selecting reac- 
tions, and a very close scrutiny of hundreds of test runs 
that were made over a period of two years. No doubt that 
some of the assumptions made had a lot of subjective intui- 
tion in them. But by and large the process of selecting the 
kinetics was made as scientifically objective as possible. 
Finally, it is very probable that even though the kinetic 
data for the four reactions in Table 5-3 appear different on 
paper, that is, each has a different activation energy and 
frequency factor, each is representative of what occurs 
within the temperature range for which the data are applica- 
ble. Therefore, the use of a reverse reaction rate constant 
tends to dampen the errors incurred in extrapolating the 
kinetic data beyond the limits set by the experimental measure 
ments. This probably contributed to the very close agreement 
in answers obtained when using different sources of data 
contrary to the author's expectations. 
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Summary of the Results for Phenolic-Nylon 

Up to this point we have analyzed the results of frozen, 
equilibrium and non-equilibrium analyses coupled to the 
virgin plastic zone. We have analyzed the energy absorbing 
mechanisms both in the virgin plastic zone and in the char 
zone . 

In the virgin plastic zone we concluded that decomposi- 
tion was the dominant energy absorbing mechanism. In the 
char zone we identified three principal energy absorbing me- 
chanisms which were transpiration cooling, chemical reactions 
and sublimation. We established that chemical reactions 
become important energy absorbing mechanisms above 2000°F, 
we pointed out that the carbon-gas and the sublimation reac- 
tions become important energy absorbing mechanisms at tempe- 
ratures above 4600°F - 4800°F. We also compared the equili- 
brium and non-equilibrium analyses boundary conditions and 
noted how each varied with increasing surface recession 
velocity. We explained why the decomposition/char zone 
temperature and mass flow rates for equilibrium and non- 
equilibrium analyses were affected by both pressures and 
surface recession velocities and how these differences affec- 
ted the shape of the heat curve for both analyses. 

We compared the three analyses at a pressure of 0.1 and 
1.0 atmospheres and at the same surface recession velocity. 

We studied the effect of both pressure and surface recession 
velocities on the total energy absorbed for both equilibrium 
and non-equilibrium. 
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We noticed that pressure had little effect on the total 
energy absorbed for equilibrium and non- equilibrium except at 
temperatures above 3600°F and 4600°F respectively. In 
analyzing the pressure effects we postulated that the carbon- 
carbon dioxide and the carbon-water reactions were probably 
the reason as to why the heat absorbed at 1.0 atmosphere was 
initially greater than at 0.1 atmosphere. We arrived at 
this conclusion by analyzing the plots of species composition 
versus temperature. Finally, we explained some of the 
numerical difficulties we encountered with the non-equilibrium 
analysis because of the stiff system of equations we were 
solving. We showed a plot of CPU time versus mass flux to 
illustrate this point. 

In conclusion then we can say that frozen flow analysis 
is a very poor approximation to the total amount of energy 
absorbed in the ablator. However, equilibrium is a reasona- 
ble approximation to non-equilibrium even though it suffers 
from its inherent simplifying assumptions that the gases are 
always in equilibrium irrespective of the temperature. 

Silicone Elastomers 

To complete this research, a brief study was made of 
another ablative composite, the silicone elastomer. This 
composite was chosen because it has good ablative properties. 
Compared to the phenolic-nylon composite, the silicone elas- 
tomer is a much denser composite; 63 lb/ft^ versus 35 lb/ft^ 
for the phenolic-nylon. The char formed by the silicone 



313 


elastomers is also denser than for phenolic-nylon (29 versus 
13). Per pound of material, the nylon absorbs 1.7 times more 
heat than the silicone elastomers. This is based on computa- 
tions performed using equilibrium analysis. 

As shown in Figure 6-26 the total heat absorbed in the 
ablator is greater at 0.1 atmosphere than at 1.0 atmosphere 
which is the behavior observed for the phenolic-nylon compo- 
site also. In addition, as shown in Figure 6-26, the compa- 
rison at a surface recession velocity of 0.034 ft/sec (pv= 

2.10 lb/ft 2 -sec) for a pressure of 0.1 and 1.0 atmosphere 
shows that there is also a crossover of the two curves, with 
the 1.0 atmosphere heat curve being slightly higher than the 
non-equilibrium curve-up to a temperature of about 3100°F. 

This crossover was also observed in the phenolic-nylon case 
both for the equilibrium and non-equilibrium analyses. 
Examination of Figures 6-27 and 6-28 shows that the carbon- 
water and carbon-carbon dioxide reactions are favored at P=l. 
atmosphere. This is evident by the rapid decrease of water 
concentration and the increase in CO concentration. It also 
appears that the silicone oxygen reaction is an important one 
since silicone oxide concentration also increases and it is 
favored at the higher pressure. 

The results for silicone elastomers at two pressure 
levels and at three surface recession velocities are presented 
on Tables 6-11 and 6-12. As was the case with phenolic-nylon 
the front surface temperature at 0.1 atmosphere is lower than 
at 1. atmosphere. Also, as was the case with phenolic-nylon, 
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Figure 6-27. Equilibrium Composition of the Silicone 
Elastomer Pyrolysis Gases in the Char at 
a Pressure of 1.0 Atmospheres. 
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Figure 6-28. Equilibrium Composition of the Silicone 
Elastomer Pyrolysis Gases in the Char 
at a Pressure of 0.1 Atmospheres. 




TABLE 6-11: Comparison of the Total Energy Absorbed at 

Three Surface Recession Velocities: Equilibrium 

Analysis (P=0.1 atm; Silicone Elastomer) . 
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the total energy absorbed at 0.1 atmosphere is greater than 
at 1.0 atmosphere. 

A non-equilibrium kinetic model for silicone elastomers 
was not developed. Characterization of this composite would 
have been more difficult than for phenolic-nylon. In addi- 
tion, the lack of kinetic data for silicone reactions, and 
the solid state reaction of silicone with silicone oxides 
and carbides would have made the kinetic analysis unreliable. 
Furthermore, we already knew that the equilibrium and non- 
equilibrium analyses for the phenolic-nylon had shown very 
close agreement, and it was postulated that probably a 
similar behavior could be expected of the silicone elasto- 
mers. 

Summary 

This chapter has presented the results of the in-depth 
response of the two ablative composites, phenolic-nylon and 
silicone elastomers. The bulk of the effort was devoted to 
phenolic-nylon since it is the composite that has shown the 
better performance in high heating rate environments. For 
phenolic-nylon we analyzed, based on the data available, 
the important energy absorbing mechanisms both for the virgin 
plastic and for the char zone. We also compared the results 
of frozen-equilibrium and non-equilibrium graphically, at two 
pressure levels, and at two surface recession velocities. 

The effects of pressure and of surface recession veloci- 
ties on the total rate of heat absorbed were studied both 
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for the phenolic-nylon and silicone elastomer composites. 

One fact should be underlined, and that is that below 3000° 
F, chemical reactions are not the most important energy ab- 
sorbing mechanisms. This is the reason why April (1) 
concluded that equilibrium analysis was not representative 
of what occurred in an ablator with a front surface tempera- 
ture of below 3000°F. However, in our studies which have 
been carried out to temperatures equal to or greater than 
5500°F, equilibrium was found to be a reasonable approxima- 
tion to characterize the ablative composites. 
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CHAPTER VII 


SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 


Summary 

This research has delat with the analysis of the heat 
absorbed of phenolic-nylon and silicone elastomer ablative 
composites. The bulk of the research was concentrated on 
phenolic-nylon since it has shown the better performance in 
high heating rate environments. Frozen, equilibrium and non- 
equilibrium analysis calculations were performed for this 
composite. The equations of change applicable to each of 
these analyses were described in Chapter III and the exten- 
sive thermodynamic and physical property data required to 
solve these equations are reported in detail in the 
Appendices. In addition, a thorough analysis and screenings 
of hundreds of reactions were done to arrive at the most 
representative set of kinetic equations for the gases re- 
sulting from the decomposition of the phenolic-nylon resin 
composite. One special precaution taken in this research 
was to calculate the equilibrium constant for each of the 
kinetic equations selected to represent the reacting pyro- 
lysis gases in the char. This was done because we wanted 
to assure ourselves that in extrapolating the kinetic 
data to temperatures above these for which the data were 
collected we did not violate any equilibrium consideration. 

The basic characteristic of this research is 
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the coupling of the polymer type virgin plastic 
zone to the char zone analysis. This was performed by 
taking into account the decomposition of the plastic compo- 
sites into gases, with their corresponding heat absorption 
and the transpiration cooling and chemical reaction effects 
that these gases have on the heat absorption in the char. 

The silicone elastomer study was performed only using 
chemical equilibrium analysis because, as explained in the 
previous chapter, essentially no kinetic data existed for 
the possible silicone reactions. The equilibrium analysis of 
the silicone elastomers demonstrates the generality of the 
program. Therefore, given the necessary input information 
as explained in Appendices A and B, any char forming ablator 
can be analyzed using the program developed during this 
research. 

Conclusions 

Based on the results of this research, the following 
conclusions are drawn: 

1. The reacting flow of pyrolysis products from 
a 40 percent nylon, 60 percent phenolic-resin 
composite is accurately described by a non- 
equilibrium model employing 15 reactions and 
19 chemical species. 

2. The total heat absorbed as predicted by the 
equilibrium and non-equilibrium analyses is a 
strong function of surface recession velocity 
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but a weak function of pressure, especially 
at temperatures below 3600°F for equilibrium 
and 4600°F for non-equilibrium. 

3. Surface recession velocity (or total mass flux) , 
which is a way of simulating a high heating 
environment, has a large effect in the decom- 
position/char zone temperature interface and in 
the total heat absorbed. 

4. Differences in the conditions at the decomposi- 
tion/char zone temperature interface between 
equilibrium and non-equilibrium are due to the 
manner in which the gas mass flux is matched 
for equilibrium. 

5. The non-equilibrium analysis for a given sur- 
face recession velocity predicts, always, a 
higher mass flux at the decomposition/ char 
zone interface because the plastic composite is 
allowed to degrade to the experimentally ' 
determined density of the char. This density is 
lower than that required to match the conditions 
for the equilibrium analysis; 13 versus r>^ 21- 
22 lbs/ft 3 . 

6 . Because the decomposition/char zone temperature 
interface is always higher than 2000°F (for all 
cases analyzed) , the gases that enter react very 
quickly. This is evident by the sharp slopes 
shown for non-equilibrium heat curves in Chapter 
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VI. 

7. As a consequence of the sharp reaction rates 
observed in the non-equilibrium (or kinetic's) 
case, the need to reduce the Runge-Kutta step 
size was necessary to maintain the stability 
of the numerical solution. This resulted in 
excessive amounts of computer time which limi- 
ted the number of cases analyzed. The lowest 
surface recession velocity analyzed was for 
v=0.02 ft/sec ?v=0.35 lb/ft^-sec) . No such 
restrictions were necessary for equilibrium or 
frozen. 

8. The principal energy absorbing mechanism in the 
plastic region is the decomposition process. 

9. In the char the three principal energy absorbing 
mechanisms are transpiration cooling, heat ab- 
sorbed by the reacting gases, and that due to 
sublimation. 

10. The equilibrium and non-equilibrium heat curves 
diverge at the higher temperature because of 
carbon sublimation. 

11. The numerical difficulties experienced with the 
non-equilibrium analysis was due to a phenomenon 
called stiffness which is caused by very fast 
reactions. These fast reactions usually occur 
at temperatures above 3000°F. 

12. The use of reverse reaction rate constants in the 
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non-equilibrium analysis prevented the 
violation of equilibrium constraints that 
might have occurred in extrapolating the 
kinetic data, beyond the limit set by the 
experimental measurements . 

In-Depth Analysis of Silicone Elastomers ; The mathe- 
matical formulation for the in-depth analysis of silicone 
elastomers was the same as for phenolic-nylon, with the 
exception that the degradation phenomena of the virgin or 
plastic zone was described by a single degradation reaction. 
The total heat absorbed by this composite was calculated 
using equilibrium analysis. Two pressure levels, and mass 

fluxes similar to those used for phenolic -nylon (0.7, 1.05, 

2 

etc., lbs/ft -sec) were used. As we saw xn Chapter VI the 
total amount of heat absorbed below 3000°F showed the same 
insensitivity to pressure as phenolic-nylon. Similar to the 
case of phenolic-nylon, surface recession velocity has a 
marked effect in the total heat absorbed by the ablator. 

In conclusion, silicone elastomers show the same 
qualitative trends observed in phenolic-nylon for the several 
parameter studies considered. However, phenolic-nylon is a 
more efficient ablator because per pound of material it 
absorbs more heat. 

Recommendations 


The following are general recommendations based on the 
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results of our research. 

Thermodynamics : In the area of thermodynamics, the free 

energy and heat capacity data for over 90 components was 
compared (where this comparison was possible) using data from 
NASA (1) , JANAF thermochemical tables (2) , API Project 44 

(3) , and data reported by Los Alamos Scientific Laboratory 

(4) . In most cases, the data for light gases was consistent. 
Where we had the choice, we selected the data from NASA 
because all of the components required to describe the equi- 
librium composition, except for phenol and benzene, which 
were not included in NASA's work, were in the form of an 
easily used polynomial (a+bT+cT 2 +. . ,+eT 4 ) . On the other 
hand, API Project 44 (which presented the data in tabular 
form) and that of Los Alamos, required transformation and 
data fitting to conform to NASA's form. However, for such 
high molecular weight components as benzene and phenol we 
were forced to use non-NASA data (2) . There were compounds 
such as toluene, 2, 4-xylenol that were reported by Sykes 

(5) which were not included in the thermodynamic analysis 
because of lack of thermodynamic data; although these higher 
molecular weight compounds are probably unimportant as far 

as the impact on the total energy absorbed. However, for the 
sake of thoroughness, if the data becomes available, they 
should be incorporated into the equilibrium analysis. 

Kinetics : A great deal of time and effort was devoted 

to narrow down the possible chemical reactions. In some 
cases we were fortunate to find more than one reference to 
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the final set of reactions selected. (See Table 5-2) . In 
such cases, we did conduct studies on the effect of kinetic 
data, but as shown in Table 5-3, there were only four reac- 
tions where this was possible. Therefore, a sensitivity 
study of the effect of activation energy and frequency fac- 
tor changes on the total energy absorbed should be conducted, 
especially for the phenol and benzene reactions. 

Decomposition Zone ? The analysis of the decomposition 
zone could be improved and a more accurate fit of the thermo- 
gravimetric data can be made. The fit used was preliminary 
(5) and some of it reported by private communications with 
NASA (6) . When the more accurate and final fit is available 
it should be used. These comments apply also to the silicone 
elastomers since the data used was based on some preliminary 
results . 

Runge-Kutta : The method used to solve the equations of 

change was a standard fourth order Runge-Kutta. This method 

worked well for non-equilibrium analysis for temperatures 

o . . 

below 3000 F. However, above this temperature numerical 

difficulties were encountered because of the very rapid 
reactions which occurred; i.e., stiffness of the set of 
equations. To be able to get around this stiffness of the 
set of equations, it was necessary to reduce the step size 
by four orders of magnitude from that of the equilibrium ana- 
lysis. 

Therefore, in future studies of this kind, implicit 
techniques should be used to avoid the excessive use of com- 



329 


puter time and avoid the numerical difficulties encountered 
in this research. 

Free Energy Minimization Formulation ; In the carbon- 
hydrogen-oxygen-nitrogen system (that of the phenolic-nylon) 
the matrix of coefficient required to solve the equilibrium 
of an all gaseous mixture is always five (see Chapter IV for 
details) . However# if solids are to be considered in the 
formulation# the rank of the matrix of coefficients is in- 
creased by one every time a solid is added to the formulation. 
For phenolic -nylon the rank of the matrix was increased by 
one since carbon was the only solid species considered. 
However, the rank of the matrix was increased by two for the 
silicone elastomers since both carbon and silicone were the 
solid species. 

The formulation just described creates a problem in that 
when the concentration of the solid species decreases or 
vanishes (as when carbon begins to sublime and the solid car- 
bon disappears) , the matrix becomes singular. The singula- 
rity is due to the fact that all the coefficients of the 
row-column combination of the solid species approach zero. 

For the phenolic-nylon case# the computations were stopped 
when the concentration of carbon was less than 10” 2 since it 
was considered that at this concentration no more char exis- 
ted. This of course, did not result in any serious error. 

However, in the silicone elastomer case the same criteria 
was used; i.e., when any of the solids reached a concentra- 
tion of 10 -2 # the temperature at which this occurred was 
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assumed to be the front end of the char. Unfortunately, 
some carbon still remained (about 1/2 mole of carbon per 
mole of gas) in the silicone analysis. Therefore, an automa- 
tic procedure should be implemented to reduce the rank of the 
matrix by one every time a solid reaches a very small concen- 
tration so as to avoid the problem encountered with the 
silicone case. 

Sublimation : In future work, a sublimation model to 

account for the carbon solid to carbon gas phase change 
should be incorporated into the non-equilibrium analysis 
because, although small, there was a residual amount of 
solid carbon that remained at the front end of the char which 
should completely disappear. 
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NOMENCLATURE 


SYMBOL 

DESCRIPTION 

UNITS 

A 

area 

L 2 

A. 

l 

frequency factor 


A. . 
ID 

identification of species in 
the jth chemical reaction. 

dimensionless 

a 

ration of char density to vir- 
gin material density (see 
Equation (2-23) . 

dimensionless 

a . . 
iD 

gram atoms of element j per 
mole of species i. 

atoms/mole 

B . 
X 

Runge-Kutta parameters in the 
species continuity numerical 
solution method. 

M/L 2 t 

b. 

D 

moles of element j defined by 
Equation (4-1) . 

moles 

b' . 
D 

moles of element j in the gas 
phase (Equation (4-34) . 

moles 

C . 

Prl 

heat capacity of component 
i at constant pressure. 

L 2 /t 2 T 

c 

P'9 

average heat capacity of a 
gas mixture . 

L 2 /t 2 T 

C k 

concentration of component K. 

3 

moles/L 

c 

number of composites in the 
virgin material. 

dimensionless 

c . 

1 

free energy function defined for 
species i defined by Equation 
(4-12). 

dimensionless 

E 

energy of activation. 

ML 2 /mole-t 2 



DESCRIPTION 


UNITS 


non-dimensional energy of 
activation. 

mathematical function defined 
by Equation (3-39). 

free energy function defined 
by Equation (4-13). 

molar free energy. 

molar free energy at standard 
state of 298 K and 1 atm. 

free energy function defined 
by Equation (4-11). 

any mathematical function. 

fugacity of species i. 

fugacity at standard state. 


dimensionless 

T/L 2 

ML 2 /t 2 

ML 2 /t 2 

ML 2 /t 2 

dimensionless 
dimensionless 
M/lt 2 ' 
M/Lt 2 


augmented function of the 
quadratic approximation to the 


free energy function defined 2 2 

by Equation (4-25). ML /t 

mathematical function defined 
by Equation (3-38). T/L 

2 

acceleration of gravity. L/t 

2 2 

enthalpy. ML /t 

integration step size. L 

2 

molar flux. moles/L t 

2 

mass flux. M/L t 

thermal conductivity. ML/t^T 


t"”^(L^/moles) n_ ^ 


reaction rate constant 



DESCRIPTION 


effective thermal conduc- 
tivity. 

ablator thickness. 

total number of chemical 
reactions. 

total number of chemical 
species, (gases+solids) . 

molecular weight. 


total number of chemical 
elements . 

molal flux of component j. 


total number of gas species, 
total pressure. 

stochiometric coefficient 
of product i in reaction 3 . 


volumetric flow rate. 


heat of pyrolysis. 

quadratic approximation 
to the free energy function 
defined by Equation (4-24). 


non-dimensional heat of pyro- 
lysis defined by Equation 
(2-40). 


energy transfer by conduction, 
convection or radiation. 


UNITS 

3 

ML/t T 
L 

dimensionless 

dimensionless 

M/moles 

dimensionless 

2 

moles/L t 

dimensionless 

M/Lt 2 

dimensionless 

L 3 /t 

M/Lt 3 

ML 2 /t 2 

dimensionless 

M/Lt 3 
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SYMBOL DESCRIPTION 


energy absorbed by the 
degradation of plastic com- 
posites, Equation (3-19). 


R 


ideal gas constant. 


chemical reaction rate. 


K effective chemical reaction 

rate for gas and solid species 
defined in Equation (3-22). 



stochiometric coefficient 
of reactant i in reaction j. 


S power on the temperature in 

the functional expression to 
calculate reaction rate cons- 
tants in Equation (5-3) 


number of solid species. 

T temperature • 

T temperature gradient. 

T - characteristic time as defined 

by Equation (5-9). 

T relaxation time as defined 

by Equation (5-10). 


t 


time. 


UNITS 

M/Lt 5 

Ml 2 /t 2 T-mole 

mole/l^t 

mole/l t 

dimensionless 

dimensionless 

dimensionless 

T 

T/L 

t 

t 

t 


u 


gas velocity 


L/T 
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SYMBOL DESCRIPTION UNITS 

3 

v volume . L 


v 


surface recession velocity. L/t 


W 


total mass flux. 


M/L 2 t 


W mass flux based on the 

p velocity in the pore’s 2 

space. M/1 t 


X. conversion of species j as 

** defined by Equation (5-6). dimensionless 


calculated value of moles of 
species i in the free energy 
minimization calculations. moles 


x total moles of gaseous moles 

species. 

£ distance defined by Equation 

(2-29). I 

assumed value of moles of 
species i in the free energy 
minimization calculations. mole 

a 

Y distance defined by Equation 

(2-6). L 

z distance in the axial 

direction. 1 
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SYMBOL 

z 

Greek : 
a 


6 


Y 


6 

A 

V 

e 

e 

n 

9 

X 

y 

n 

p 


DESCRIPTION 

UNITS 

compressibility factor. 

dimensionless 

viscous coefficient in the 
modified Darcy's Law Equation. 

L- 2 

constant defined by Equation 
(2-24) 

dimensionless 

activity coefficient defined 
Equation (4-4) . 

dimensionless 

inertial coefficient in the 
modified Darcy's Law Equation. 

L" 1 

non-dimensional energy of acti- 
vation defined by Equation 
(2-45). 

dimensionless 

permiability of a porous medium. 

x . 2 

non-dimensional heat capacity 
defined by Equation (2-39) . 

dimensionless 

Kronocker Delta. 

dimensionless 

a difference between two 
parameters . 

dimensionless 

del operator. 

dimensionless 

porosity. 

dimensionless 

emissivity 

dimensionless 

dimensionless char 
distance . 

dimensionless 

dimensionless temperature. 

dimensionless 

eigen value. 

dimensionless 

viscosity. 

M/Lt 

Lagrange multiplier. 

ML 2 /t 2 

density. 

M/L 3 


SYMBOL 


DESCRIPTION 


UNITS 


a 

a 

T 

Subscript: 

c 

e 

g 

L 

o 

P 

r 

s 

T 


Stefan-Boltzman constant, 
collision diameter, 
shear stress. 

convection or conduction. 

effective . 

gas. 

front surface of char. 

initial condition. 

pressure. 

pyrolysis . 

radiation. 

solid. 

temperature . 


Superscript: 

o standard or reference state. 


M/t 3 T 

L 

M/Lt 2 


derivative 



